!  LES2LWC
! Converts Franks LWC netcdf output file from the UCLA LES with 3D LWC
! fields for several times to several of Frank's ascii LWC fields.  
! Only outputs part of the height range with max LWC above a specified 
! threshold to the ascii LWC fields.  
!   Also outputs the LES field grid size and spacing, the cloud fraction 
! and mean cloudy column liquid water path, and profiles of cloud fraction
! and mean cloud LWC.
!
!
!  g95 -O2 -o les2lwc les2lwc.f90  -lnetcdf 
!      -I/usr/local/netcdf-3.6.3/include/  -L/usr/local/netcdf-3.6.3/lib/
!
!   Frank Evans   July 2009

MODULE LES2LWC_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, 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)
  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, '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 CALCULATE_CLOUDFRAC_LWP (NX, NY, NZ, NT, LWPTHRESH, LWCTHRESH, &
                                    ZLEVELS, LWC, IZMIN, IZMAX, &
                                    CLDFRACPROF, CLDLWCPROF)
 ! Calculates and prints out the cloud fraction and liquid water path,
 ! and the height range that contains the LWC levels above the
 ! LWC threshold.  The input LWP threshold is used to define the vertical
 ! cloud fraction, and the mean cloud LWP is calculated.  The input LWC
 ! threshold is use to define the cloud fraction for each level, and
 ! the profiles of cloud fraction and mean cloud LWC are calculated.
 ! The max LWC profile of the cloudy columns is computed, and then 
 ! searched for the cloudy height range using the input LWC threshold.  
 ! The indices of this height range are returned for each time.
  IMPLICIT NONE
  INTEGER, INTENT(IN)  :: NX, NY, NZ, NT
  REAL,    INTENT(IN)  :: LWPTHRESH, LWCTHRESH, ZLEVELS(NZ)
  REAL,    INTENT(IN)  :: LWC(NZ,NX,NY,NT)
  INTEGER, INTENT(OUT) :: IZMIN(NT), IZMAX(NT)
  REAL,    INTENT(OUT) :: CLDFRACPROF(NZ,NT), CLDLWCPROF(NZ,NT)
  INTEGER :: J, IZ, NCLD
  REAL    :: CLDFRAC, AVGCLDLWP
  REAL, ALLOCATABLE :: LWP(:,:), LWCPROF(:)

  ALLOCATE (LWP(NX,NY), LWCPROF(NZ))
  PRINT '(A)', 'Itime CldFrac AvgCldLWP  Zmin  Zmax'
  DO J = 1, NT
     ! Compute the 2D LWP field
    LWP(:,:) = 0.0
    DO IZ = 1, NZ-1
      LWP(:,:) = LWP(:,:) + 0.5*(LWC(IZ,:,:,J)+LWC(IZ+1,:,:,J)) &
                    *1000.0*(ZLEVELS(IZ+1)-ZLEVELS(IZ))
    ENDDO
     ! Use the LWP threshold to calculate the cloud fraction and average LWP
    NCLD = COUNT(LWP(:,:)>LWPTHRESH)
    CLDFRAC = FLOAT(NCLD)/(NX*NY)
    IF (NCLD > 0) THEN
      AVGCLDLWP = SUM(LWP(:,:), LWP(:,:)>LWPTHRESH)/NCLD
    ELSE
      AVGCLDLWP = 0.0
    ENDIF

     ! Compute the profiles of cloud fraction and mean cloud LWC
    DO IZ = 1, NZ
      NCLD = COUNT(LWC(IZ,:,:,J)>LWCTHRESH)
      CLDFRACPROF(IZ,J) = FLOAT(NCLD)/(NX*NY)
      IF (NCLD > 0) THEN
        CLDLWCPROF(IZ,J) = SUM(LWC(IZ,:,:,J), LWC(IZ,:,:,J)>LWCTHRESH) /NCLD
      ELSE
        CLDLWCPROF(IZ,J) = 0.0
      ENDIF
    ENDDO

     ! Calculate the max LWC profile
    DO IZ = 1, NZ
      LWCPROF(IZ) = MAXVAL(LWC(IZ,:,:,J))
    ENDDO
     ! Find the height range that includes the max LWC above
     !   the specified threshold
    IZMAX(J) = 1
    IZMIN(J) = 1
    IZ = 1
    DO WHILE (IZ<=NZ)
      IF (LWCPROF(IZ) > LWCTHRESH) IZMAX(J) = IZ
      IZ = IZ + 1
    ENDDO
    IZ = NZ
    DO WHILE (IZ>=1)
      IF (LWCPROF(IZ) > LWCTHRESH) IZMIN(J) = IZ
      IZ = IZ - 1
    ENDDO
     ! Print out the results for this time
    PRINT '(2X,I2,3X,F5.3,2X,F6.2,4X,2(1X,F5.3))', &
           J, CLDFRAC, AVGCLDLWP, ZLEVELS(IZMIN(J)), ZLEVELS(IZMAX(J))
  ENDDO
  DEALLOCATE (LWP, LWCPROF)
END SUBROUTINE CALCULATE_CLOUDFRAC_LWP



SUBROUTINE WRITE_PROFILE_FILE (PROFFILE, NZ, NT, INCDFFILE, LWCTHRESH, &
                               ZLEVELS, CLDFRACPROF, CLDLWCPROF)
                               
  IMPLICIT NONE
  CHARACTER(LEN=*), INTENT(IN)  :: PROFFILE, INCDFFILE
  INTEGER, INTENT(IN)           :: NZ, NT
  REAL,    INTENT(IN)           :: LWCTHRESH, ZLEVELS(NZ)
  REAL,    INTENT(IN)           :: CLDFRACPROF(NZ,NT), CLDLWCPROF(NZ,NT)
  INTEGER :: J, K
  LOGICAL :: CLOUDY(NT)

  CLOUDY(:) = (SUM(CLDFRACPROF(:,:),DIM=1) > 0.0)
  OPEN (UNIT=1, FILE=PROFFILE, STATUS='UNKNOWN')
  WRITE (1, '(A)') '! Cloud fraction and mean cloud liquid water content profiles.'
  WRITE (1, '(A,A60)') '! LES LES netcdf file: ',INCDFFILE
  WRITE (1, '(A,F7.5)') '! LWC threshold (g/m^3): ', LWCTHRESH
  WRITE (1, '(A)') '! Height CldFrac AvgLWC  CldFrac AvgLWC ...'
  WRITE (1,'(A8)',ADVANCE='NO') '! Times:'
  DO J = 1, NT
    IF (CLOUDY(J)) THEN
      WRITE (1,'(5X,I3,5X)',ADVANCE='NO') J
    ENDIF
  ENDDO
  WRITE (1,*)
  DO K = 1, NZ
    WRITE (1,'(1X,F7.4)',ADVANCE='NO') ZLEVELS(K)
    DO J = 1, NT
      IF (CLOUDY(J)) THEN
        WRITE (1,'(2X,F5.3,1X,F5.3)',ADVANCE='NO') CLDFRACPROF(K,J),CLDLWCPROF(K,J)
      ENDIF
    ENDDO
    WRITE (1,*)
  ENDDO
  CLOSE (1)
END SUBROUTINE WRITE_PROFILE_FILE



SUBROUTINE WRITE_LWC_FILES (OUTFILEBASE, NX, NY, NZ, NT, IZMIN, IZMAX, & 
                           DROPCONC, DELX, DELY, ZLEVELS, TEMPS, LWC)
 ! Outputs the two parameter liquid water content files for all the NT
 ! times in the LWC field.  Only outputs the height levels 
 ! specified by the IZMIN and IZMAX index ranges.  Only outputs those
 ! grid cells with LWC above MINLWC.  
  IMPLICIT NONE
  CHARACTER(LEN=*), INTENT(IN)  :: OUTFILEBASE
  INTEGER, INTENT(IN)           :: NX, NY, NZ, NT
  INTEGER, INTENT(INOUT)        :: IZMIN(NT), IZMAX(NT)
  REAL,    INTENT(IN)           :: DELX, DELY, ZLEVELS(NZ), TEMPS(NZ,NT)
  REAL,    INTENT(IN)           :: LWC(NZ,NX,NY,NT), DROPCONC
  REAL, PARAMETER :: MINLWC=0.0001
  INTEGER :: J, K, IX, IY, IZ, NZC
  REAL, PARAMETER :: ALPHA=7.0 ! gamma distribution width parameter
  REAL    :: REFF
  CHARACTER(LEN=80) :: FILENAME

  DO J = 1, NT
    IF (IZMIN(J) == 1 .AND. IZMAX(J) == 1) THEN
      IZMIN(J) = NZ/3 ; IZMAX(J) = NZ/3
    ENDIF
    NZC = IZMAX(J) - IZMIN(J) + 1
    IF (NZC > 0) THEN
      WRITE (FILENAME, '(A,I2.2,A)') TRIM(OUTFILEBASE), J, '.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,'(100(1X,F7.4))') ZLEVELS(IZMIN(J):IZMAX(J))
      WRITE (1,'(100(1X,F6.2))') TEMPS(IZMIN(J):IZMAX(J),J)
      DO IX = 1, NX
        DO IY = 1, NY
          DO IZ = 1, NZC
            K = IZ + IZMIN(J) - 1
            IF (LWC(K,IX,IY,J) > MINLWC) THEN
              REFF = 100* ( LWC(K,IX,IY,J) *(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,J), REFF
            ENDIF
          ENDDO
        ENDDO
      ENDDO
      CLOSE (1)
    ENDIF
  ENDDO
END SUBROUTINE WRITE_LWC_FILES




SUBROUTINE WRITE_ATMOS_PROF_FILES (OUTFILEBASE, NX, NY, NZ, NT, & 
                                   ZLEVELS, PRES, TEMPS, VAPMIX)
 ! Outputs the atmospheric profile files for all the NT times in the 
 ! LWC field.  
 !  Format: height (km), pressure (mb), temperature (K), water vapor (g/g)
  IMPLICIT NONE
  CHARACTER(LEN=*), INTENT(IN)  :: OUTFILEBASE
  INTEGER, INTENT(IN)           :: NX, NY, NZ, NT
  REAL,    INTENT(IN)           :: ZLEVELS(NZ)
  REAL,    INTENT(IN)           :: PRES(NZ,NT), TEMPS(NZ,NT), VAPMIX(NZ,NT)
  INTEGER :: J, K, IZ
  REAL    :: F
  CHARACTER(LEN=80) :: FILENAME

  DO J = 1, NT
    WRITE (FILENAME, '(A,I2.2,A)') TRIM(OUTFILEBASE), J, '.atmos'
    OPEN (UNIT=2, FILE=FILENAME, STATUS='UNKNOWN')
    F = (0.0-ZLEVELS(1))/(ZLEVELS(2)-ZLEVELS(1))
    WRITE (2,'(F7.4,1X,F7.2,1X,F6.2,1X,E10.4)') &
       0.0, (1-F)*PRES(1,J)+F*PRES(2,J), (1-F)*TEMPS(1,J)+F*TEMPS(2,J), &
       0.001*((1-F)*VAPMIX(1,J)+F*VAPMIX(2,J))
    DO IZ = 2, NZ
       WRITE (2,'(F7.4,1X,F7.2,1X,F6.2,1X,E10.4)') &
         ZLEVELS(IZ), PRES(IZ,J), TEMPS(IZ,J), 0.001*VAPMIX(IZ,J)
    ENDDO
    CLOSE (2)
  ENDDO
END SUBROUTINE WRITE_ATMOS_PROF_FILES

END MODULE LES2LWC_SUB



PROGRAM LES2LWC
  USE LES2LWC_SUB
  IMPLICIT NONE
  REAL                 :: LWPTHRESH, LWCTHRESH, DROPCONC
  CHARACTER(LEN=72)    :: INCDFFILE, OUTFILEBASE, PROFFILE
  INTEGER              :: NX, NY, NZ, NT
  LOGICAL              :: ATMOSOUTFLAG
  REAL                 :: DELX, DELY
  INTEGER, ALLOCATABLE :: IZMIN(:), IZMAX(:)
  REAL,    ALLOCATABLE :: ZLEVELS(:), PRES(:,:), TEMPS(:,:), VAPMIX(:,:)
  REAL,    ALLOCATABLE :: LWC(:,:,:,:)
  REAL,    ALLOCATABLE :: CLDFRACPROF(:,:), CLDLWCPROF(:,:)

  INTEGER :: I, K

  WRITE (*,*) 'Input UCLA LES LWC netcdf file name'
  READ (*,'(A)') INCDFFILE
  WRITE (*,*) INCDFFILE

  WRITE (*,*) 'Output cloud profile file (or NONE)'
  READ (*,'(A)') PROFFILE
  WRITE (*,*) PROFFILE

  WRITE (*,*) 'Base name for output ascii LWC and atmos profile files (or NONE)'
  READ (*,'(A)') OUTFILEBASE
  WRITE (*,*) OUTFILEBASE

  WRITE (*,*) 'Output atmospheric profile files (T or F)'
  READ (*,*) ATMOSOUTFLAG
  WRITE (*,*) ATMOSOUTFLAG

  WRITE (*,*) 'LWP threshold for cloud fraction (g/m^2)'
  READ (*,*) LWPTHRESH
  WRITE (*,*) LWPTHRESH

  WRITE (*,*) 'LWC (g/m^3) threshold for outputing layers'
  READ (*,*) LWCTHRESH
  WRITE (*,*) LWCTHRESH

  WRITE (*,*) 'Cloud droplet number concentration (cm^-3)'
  READ (*,*) DROPCONC
  WRITE (*,*) DROPCONC


  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))
  ALLOCATE (IZMIN(NT), IZMAX(NT))
  ALLOCATE (CLDFRACPROF(NZ,NT), CLDLWCPROF(NZ,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)

  CALL READ_CDF_LWC (INCDFFILE, NX, NY, NZ, NT, LWC, PRES, TEMPS, VAPMIX)

  CALL CALCULATE_CLOUDFRAC_LWP (NX, NY, NZ, NT, LWPTHRESH, LWCTHRESH, &
                           ZLEVELS, LWC, IZMIN, IZMAX, CLDFRACPROF, CLDLWCPROF)

  IF (PROFFILE /= 'NONE') THEN
    CALL WRITE_PROFILE_FILE (PROFFILE, NZ, NT, INCDFFILE, LWCTHRESH, &
                             ZLEVELS, CLDFRACPROF, CLDLWCPROF)
  ENDIF

  IF (OUTFILEBASE /= 'NONE') THEN
    CALL WRITE_LWC_FILES (OUTFILEBASE, NX, NY, NZ, NT, IZMIN, IZMAX, & 
                         DROPCONC, DELX, DELY, ZLEVELS, TEMPS, LWC)
  ENDIF

  IF (ATMOSOUTFLAG) THEN
    CALL WRITE_ATMOS_PROF_FILES (OUTFILEBASE, NX, NY, NZ, NT, & 
                                 ZLEVELS, PRES, TEMPS, VAPMIX)
  ENDIF

  DEALLOCATE (ZLEVELS, PRES, TEMPS, VAPMIX, LWC)
  DEALLOCATE (IZMIN, IZMAX, CLDFRACPROF, CLDLWCPROF)
END PROGRAM LES2LWC

