! LES2LWC2 ! Converts Franks new LWC netcdf output file from the UCLA LES with 3D LWC ! and relative humidity fields for several times to several of Frank's ! ascii cloud LWC fields (old type 2) with either LWC and derived effective ! radius or LWC and relative humidity. ! ! gfortran -O2 -o les2lwc2 les2lwc2.f90 -lnetcdf -lnetcdff ! -I/usr/local/include/ -L/usr/local/lib/ ! ! Frank Evans April 2012 MODULE LES2LWC2_SUB CONTAINS SUBROUTINE READ_CDF_SIZE (INFILE, NX, NY, NZ, NT) ! Returns the size of the LES grid (the LWC file does not include the ! four edge grid points). IMPLICIT NONE CHARACTER(LEN=*), INTENT(IN) :: INFILE INTEGER, INTENT(OUT) :: NX, NY, NZ, NT INTEGER :: status, ncid, dimid INCLUDE 'netcdf.inc' ! Get ID of netcdf file status = NF_OPEN (infile,0,ncid) if (status /= 0) stop 'Error opening CDF file' status = NF_INQ_DIMID (ncid, 'time', dimid) status = NF_INQ_DIMLEN (ncid, dimid, NT) if (status /= 0) stop 'Error getting length of time dimension' ! print *, 'nt=',NT status = NF_INQ_DIMID (ncid, 'zt', dimid) status = NF_INQ_DIMLEN (ncid, dimid, NZ) if (status /= 0) stop 'Error getting length of Z dimension' ! print *, 'nz=',NZ status = NF_INQ_DIMID (ncid, 'yt', dimid) status = NF_INQ_DIMLEN (ncid, dimid, NY) if (status /= 0) stop 'Error getting length of Y dimension' ! print *, 'ny=',NY status = NF_INQ_DIMID (ncid, 'xt', dimid) status = NF_INQ_DIMLEN (ncid, dimid, NX) if (status /= 0) stop 'Error getting length of X dimension' ! print *, 'nx=',NX status = NF_CLOSE (ncid) END SUBROUTINE READ_CDF_SIZE SUBROUTINE READ_CDF_GRID (INFILE, NX, NY, NZ, NT, DELX, DELY, ZLEVELS) ! Opens the netcdf file and reads the X, Y, Z gridpoint locations. ! Assumes constant grid spacing for X and Y and returns the spacing in ! DELX and DELY. Output units are km. IMPLICIT NONE INTEGER, INTENT(IN) :: NX, NY, NZ, NT CHARACTER(LEN=*), INTENT(IN) :: INFILE REAL, INTENT(OUT):: DELX, DELY, ZLEVELS(NZ) INTEGER :: status, ncid, arrayid REAL, ALLOCATABLE :: xvals(:), yvals(:), zvals(:) INCLUDE 'netcdf.inc' allocate (xvals(nx+4), yvals(ny+4), zvals(nz)) ! Get ID of netcdf file status = NF_OPEN (infile,0,ncid) if (status /= 0) stop 'Error opening CDF file' status = NF_INQ_VARID (ncid, 'xt', arrayid) status = NF_GET_VAR_REAL (ncid, arrayid, xvals) if (status /= 0) stop 'Error getting X locations' DELX = 0.001*(xvals(nx)-xvals(1))/(nx-1) status = NF_INQ_VARID (ncid, 'yt', arrayid) status = NF_GET_VAR_REAL (ncid, arrayid, yvals) if (status /= 0) stop 'Error getting Y locations' DELY = 0.001*(yvals(ny)-yvals(1))/(ny-1) status = NF_INQ_VARID (ncid, 'zt', arrayid) status = NF_GET_VAR_REAL (ncid, arrayid, zvals) if (status /= 0) stop 'Error getting Z levels' ZLEVELS(:) = 0.001*zvals(:) status = NF_CLOSE (ncid) deallocate (xvals, yvals, zvals) END SUBROUTINE READ_CDF_GRID SUBROUTINE READ_CDF_LWC (INFILE, NX, NY, NZ, NT, LWC, RH, PRES, TEMPS, VAPMIX) ! Opens the netcdf file and reads the liquid water content (g/m^3) field. IMPLICIT NONE INTEGER, INTENT(IN) :: NX, NY, NZ, NT CHARACTER(LEN=*), INTENT(IN) :: INFILE REAL, INTENT(OUT) :: LWC(NZ,NX,NY,NT), RH(NZ,NX,NY,NT) REAL, INTENT(OUT) :: PRES(NZ,NT), TEMPS(NZ,NT), VAPMIX(NZ,NT) INTEGER :: status, ncid, arrayid INCLUDE 'netcdf.inc' ! Get ID of netcdf file status = NF_OPEN (infile,0,ncid) if (status /= 0) stop 'Error opening CDF file' status = NF_INQ_VARID (ncid, 'lwc', arrayid) status = NF_GET_VAR_REAL (ncid, arrayid, LWC) if (status /= 0) stop 'Error getting LWC' status = NF_INQ_VARID (ncid, 'rh', arrayid) status = NF_GET_VAR_REAL (ncid, arrayid, RH) if (status /= 0) stop 'Error getting RH' status = NF_INQ_VARID (ncid, 'p1', arrayid) status = NF_GET_VAR_REAL (ncid, arrayid, PRES) if (status /= 0) stop 'Error getting PRES' status = NF_INQ_VARID (ncid, 't1', arrayid) status = NF_GET_VAR_REAL (ncid, arrayid, TEMPS) if (status /= 0) stop 'Error getting TEMPS' status = NF_INQ_VARID (ncid, 'qv1', arrayid) status = NF_GET_VAR_REAL (ncid, arrayid, VAPMIX) if (status /= 0) stop 'Error getting VAPMIX' status = NF_CLOSE (ncid) END SUBROUTINE READ_CDF_LWC SUBROUTINE WRITE_LWC_RH_FILES (OUTFILEBASE, NX, NY, NZ, NT, IT, RHOUT, & DROPCONC, DELX, DELY, ZLEVELS, TEMPS, LWC, RH) ! Outputs the two parameter liquid water content files for one (IT) of the ! NT times in the LWC field, but the effective radius can be replaced by ! the relative humidity if RHOUT=.true. IMPLICIT NONE CHARACTER(LEN=*), INTENT(IN) :: OUTFILEBASE INTEGER, INTENT(IN) :: NX, NY, NZ, NT, IT LOGICAL, INTENT(IN) :: RHOUT REAL, INTENT(IN) :: DELX, DELY, ZLEVELS(NZ), TEMPS(NZ,NT) REAL, INTENT(IN) :: LWC(NZ,NX,NY,NT), RH(NZ,NX,NY,NT), DROPCONC REAL, PARAMETER :: MINLWC=0.0001 INTEGER :: K, IX, IY, IZ, NZC REAL, PARAMETER :: ALPHA=7.0 ! gamma distribution width parameter REAL :: REFF CHARACTER(LEN=80) :: FILENAME NZC = NZ IF (NZC > 0) THEN WRITE (FILENAME, '(A,I2.2,A)') TRIM(OUTFILEBASE), IT, '.lwc' ! print '(A,A60)', 'Writing : ', FILENAME OPEN (UNIT=1, FILE=FILENAME, STATUS='UNKNOWN') WRITE (1,'(A)') '2 parameter LWC file' WRITE (1,'(3(1X,I3))') NX, NY, NZC WRITE (1,'(2(1X,F6.4))') DELX, DELY WRITE (1,'(150(1X,F7.4))') ZLEVELS(1:NZ) WRITE (1,'(150(1X,F6.2))') TEMPS(1:NZ,IT) DO IX = 1, NX DO IY = 1, NY DO IZ = 1, NZC K = IZ IF (RHOUT)THEN WRITE (1, '(3(I3,1X),F6.4,1X,F5.3)') & IX, IY, IZ, LWC(K,IX,IY,IT), RH(K,IX,IY,IT) ELSE IF (LWC(K,IX,IY,IT) > MINLWC) THEN REFF = 100* ( LWC(K,IX,IY,IT) *(0.75*(ALPHA+3)*(ALPHA+3)) & /(3.14159*(ALPHA+1)*(ALPHA+2)*DROPCONC) )**(1.0/3) WRITE (1, '(3(I3,1X),F6.4,1X,F5.2)') & IX, IY, IZ, LWC(K,IX,IY,IT), REFF ENDIF ENDIF ENDDO ENDDO ENDDO CLOSE (1) ENDIF END SUBROUTINE WRITE_LWC_RH_FILES END MODULE LES2LWC2_SUB PROGRAM LES2LWC2 USE LES2LWC2_SUB IMPLICIT NONE REAL :: DROPCONC LOGICAL :: RHOUT INTEGER :: IT1, IT2, IT CHARACTER(LEN=72) :: INCDFFILE, OUTFILEBASE INTEGER :: NX, NY, NZ, NT REAL :: DELX, DELY REAL, ALLOCATABLE :: ZLEVELS(:), PRES(:,:), TEMPS(:,:), VAPMIX(:,:) REAL, ALLOCATABLE :: LWC(:,:,:,:), RH(:,:,:,:) INTEGER :: I, K WRITE (*,*) 'Input UCLA LES LWC netcdf file name' READ (*,'(A)') INCDFFILE WRITE (*,*) INCDFFILE WRITE (*,*) 'Base name for output ascii LWC/RH files (or NONE)' READ (*,'(A)') OUTFILEBASE WRITE (*,*) OUTFILEBASE WRITE (*,*) 'Cloud droplet number concentration (cm^-3)' READ (*,*) DROPCONC WRITE (*,*) DROPCONC WRITE (*,*) 'Output relative humidity instead of effective radius (T or F) ' READ (*,*) RHOUT WRITE (*,*) RHOUT WRITE (*,*) 'Starting and ending time index in LES LWC netcdf file' READ (*,*) IT1, IT2 WRITE (*,*) IT1, IT2 CALL READ_CDF_SIZE (INCDFFILE, NX, NY, NZ, NT) ALLOCATE (ZLEVELS(NZ), PRES(NZ,NT), TEMPS(NZ,NT), VAPMIX(NZ,NT)) ALLOCATE (LWC(NZ,NX,NY,NT), RH(NZ,NX,NY,NT)) CALL READ_CDF_GRID (INCDFFILE, NX, NY, NZ, NT, DELX, DELY, ZLEVELS) PRINT * PRINT '(A,4(1X,I3),1X,F7.5,1X,F5.3)', & 'Nx, Ny, Nz, Nt, delx, Zmax: ', NX,NY,NZ,NT,DELX,ZLEVELS(NZ) IF (OUTFILEBASE /= 'NONE') THEN CALL READ_CDF_LWC (INCDFFILE, NX, NY, NZ, NT, LWC, RH, PRES, TEMPS, VAPMIX) DO IT = IT1, IT2 CALL WRITE_LWC_RH_FILES (OUTFILEBASE, NX, NY, NZ, NT, IT, RHOUT, & DROPCONC, DELX, DELY, ZLEVELS, TEMPS, LWC, RH) ENDDO ENDIF DEALLOCATE (ZLEVELS, PRES, TEMPS, VAPMIX, LWC, RH) END PROGRAM LES2LWC2