!############################################################################### !---------------------- Driver of the new microphysics ------------------------- !############################################################################### SUBROUTINE GSMDRIVE !------------------------------------------------------------------------------- !----- NOTE: Code is currently set up w/o threading! !------------------------------------------------------------------------------- !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . ! SUBPROGRAM: Grid-scale microphysical processes - condensation & precipitation ! PRGRMMR: Ferrier ORG: W/NP22 DATE: February 2001 ! 2001-04-xx Ferrier - Beta-tested version ! 2001-05-21 Ferrier - Added gradual latent heating to remove external waves ! 2001-05-30 Ferrier - Changed default to uniform maritime conditions for testing !------------------------------------------------------------------------------- ! ABSTRACT: ! * Merges original GSCOND & PRECPD subroutines. ! * Code has been substantially streamlined and restructured. ! * Exchange between water vapor & small cloud condensate is calculated using ! the original Asai (1965, J. Japan) algorithm. See also references to ! Yau and Austin (1979, JAS), Rutledge and Hobbs (1983, JAS), and Tao et al. ! (1989, MWR). This algorithm replaces the Sundqvist et al. (1989, MWR) ! parameterization. !------------------------------------------------------------------------------- ! Prior PROGRAM HISTORY LOG: ! *** Heritage as Subroutine GSCOND: ! 94-~?? ZHAO - ORIGINATOR ! 95-03-25 BLACK - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL ! 95-03-28 BLACK - ADDED EXTERNAL EDGE ! 98-11-02 BLACK - MODIFIED FOR DISTRIBUTED MEMORY ! *** Heritage as Subroutine PRECPD: ! 94-~?? ZHAO - ORIGINATOR ! 95-03-25 BLACK - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL ! 95-11-20 ABELES - PARALLEL OPTIMIZATION ! 96-03-29 BLACK - REMOVED SCRCH COMMON ! 96-07-18 ZHAO - NEW WMIN CALCULATION ! 96-09-25 BALDWIN - NEW SR CALCULATION ! 98-11-02 BLACK - MODIFICATION FOR DISTRIBUTED MEMORY !------------------------------------------------------------------------------- ! USAGE: CALL GSMDRIVE FROM MAIN PROGRAM EBU ! INPUT ARGUMENT LIST: ! NONE ! OUTPUT ARGUMENT LIST: ! NONE ! OUTPUT FILES: ! NONE ! Subprograms & Functions called: ! GSMCONST - initialize rain & ice lookup tables, read from external file; ! initialize constants ! GSMCOLUMN - cloud microphysics calculations over vertical columns ! UNIQUE: NONE ! LIBRARY: NONE !--- COMMON BLOCKS (input for microphysics): ! CTLBLK, LOOPS, MASKS, PHYS, VRBLS, CLDWTR, PVRBLS, ACMCLH, PPTASM, C_FRACN !--- COMMON BLOCKS ("triggers" for microphysics & statistics): ! CMICRO_START, CMICRO_STATS ! ATTRIBUTES: ! LANGUAGE: FORTRAN 90 ! MACHINE : IBM SP !------------------------------------------------------------------------ INCLUDE "parmeta.f90" INCLUDE "parm.tbl.f90" INCLUDE "mpp.h" INCLUDE "mpif.h" #include "sp.h" !---------------------------------------------------------------------- INTEGER, PARAMETER :: JAM=6+2*(JM-10), LP1=LM+1 !----------------------------------------------------------------------- LOGICAL :: RUN,FIRST,RESTRT,SIGMA,NOZ !---------------------------------------------------------------------- INCLUDE "COMM_CTLBLK.f90" INCLUDE "COMM_LOOPS.f90" INCLUDE "COMM_MASKS.f90" INCLUDE "COMM_PHYS.f90" INCLUDE "COMM_VRBLS.f90" INCLUDE "COMM_CLDWTR.f90" INCLUDE "COMM_PVRBLS.f90" INCLUDE "COMM_ACMCLH.f90" INCLUDE "COMM_PPTASM.f90" INCLUDE "COMM_C_FRACN.f90" !------------------------------------------------------------------------- !----- Key parameters passed to column microphysics (COLUMN_MICRO) ------ !------------------------------------------------------------------------- !--- Flag from INIT.F at start of model run, used in initiating statistics COMMON /CMICRO_START/ MICRO_START LOGICAL :: MICRO_START !--- This variable is for debugging purposes (if .true.) LOGICAL, PARAMETER :: PRINT_diag= .TRUE. !--- The following variables are for microphysical statistics (non-essential) INTEGER, PARAMETER :: ITLO=-60, ITHI=40, ITHILO=ITHI-ITLO+1, & ITHILO_N=ITHILO*4, ITHILO_QM=ITHILO*5, ITHILO_QT=ITHILO*22 COMMON /CMICRO_STATS/ NSTATS(ITLO:ITHI,4), QMAX(ITLO:ITHI,5), & QTOT(ITLO:ITHI,22) INTEGER :: NSTATS, NSTATS_0(ITLO:ITHI,4) REAL :: QMAX, QTOT, QMAX_0(ITLO:ITHI,5),QTOT_0(ITLO:ITHI,22) REAL, SAVE :: Thour_print, & PRECmax(2),PRECtot(2),PRECmax_0(2),PRECtot_0(2) REAL, PARAMETER :: DThour_print=3. ! Print statistics every 3 h ! REAL, PARAMETER :: DThour_print=1. ! Print statistics every 1 h ! REAL, PARAMETER :: DThour_print=0. ! Print statistics every time step !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~~~~ BEGIN section on hydrometeor fractions !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~~~~ Saved values use REAL (REAL*4) arrays rather than INTEGER*2 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ REAL :: F_ice, F_rain, F_RimeF, Fice, Frain, DUM !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~~~~ Saved values use INTEGER*2 arrays rather than REAL*4 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !--- Parameters used in converting from real to integer*2 fields for ! fraction of ice (F_ice), fraction of rain (F_rain), and mass ! ratio of rimed ice ("rime factor", F_RimeF). !! REAL, PARAMETER :: FSCALE=10000., RFSCALE=1./FSCALE, !! & RIMESCALE=100., RIMEROUND=.5/RIMESCALE !! INTEGER*2, PARAMETER :: IFMIN=0, IFMAX=FSCALE, IRimeF=RIMESCALE !! !! INTEGER*2 :: F_ice, F_rain, F_RimeF, Fice, Frain, I2DUM !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~~~~ END section on hydrometeor fractions !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !----------------------------------------------------------------------- !-------------- Local arrays & parameters in GSMDRIVE ----------------- !----------------------------------------------------------------------- !---- Comments on 23 August 2001 ! * EPSQ=1.E-20 is the lower limit for specific humidity and cloud ! condensate. The value of EPSQ will need to be changed in the other ! subroutines in order to make it consistent throughout the Eta code. REAL, PARAMETER :: EPSQ=1.E-20, GRAV=9.81, RHOL=1000., T0C=273.15, & T_ICE=-10., T_ICEK=T0C+T_ICE, RRHOL=1./RHOL, CLIMIT=1.E-12 REAL :: ARAIN, ASNOW, P_col(LM), QI_col(LM), QR_col(LM), & QV_col(LM), QW_col(LM), RimeF_col(LM), T_col(LM), THICK_col(LM), & WC_col(LM) !------------------------------------------------------------------------ !####################################################################### !########################## Begin Execution ############################ !####################################################################### !------------------------------------------------------------------------ !---------------------- Microphysical constants ------------------------- !------------------------------------------------------------------------ DTPH=DTQ2 ! Physics time step (s) !------------------------------------------------------------------------ !--------------- Initialize constants for statistics -------------------- !------------------------------------------------------------------------ IF (MICRO_START) THEN MICRO_START= .FALSE. ! No need to calculate these parameters again !--- Begin: Initialize arrays if not initialized in restart files (1/16/02) DUM=0. DO J=MYJS2,MYJE2 DO I=MYIS,MYIE IF (HBM2(I,J) > 0.5) THEN LSFC=LMH(I,J) ! "L" of surface DO L=1,LSFC DUM=MAX(DUM,F_rain(L,I,J),F_ice(L,I,J),F_RimeF(L,I,J)) IF (DUM > 0.) GO TO 110 ENDDO ENDIF ENDDO ENDDO 110 IF (DUM <= 0.) THEN IF (MYPE == 0) & WRITE(0, "(a,2I3)" ) 'Initialize internal microphysics arrays' DO J=MYJS2,MYJE2 DO I=MYIS,MYIE IF (HBM2(I,J) > 0.5) THEN LSFC=LMH(I,J) ! "L" of surface DO L=1,LSFC F_rain(L,I,J)=0. IF (T(I,J,L) <= T_ICEK) THEN F_ice(L,I,J)=1. ELSE F_ice(L,I,J)=0. ENDIF F_RimeF(L,I,J)=1. ENDDO ENDIF ENDDO ENDDO ENDIF CALL GSMCONST (DTPH) ! Initialize lookup tables & constants Thour_print=-DTPH/3600. IF (PRINT_diag) THEN !-------- Total and maximum quantities NSTATS=0 ! Microphysical statistics dealing w/ grid-point counts QMAX=0. ! Microphysical statistics dealing w/ hydrometeor mass QTOT=0. ! Microphysical statistics dealing w/ hydrometeor mass PRECmax=0. ! Maximum precip rates (rain, snow) at surface (mm/h) PRECtot=0. ! Total precipitation (rain, snow) accumulation at surface ENDIF ENDIF !----------------------------------------------------------------------- !------------------- Loop horizontally through domain ------------------ !----------------------------------------------------------------------- DO 100 J=MYJS2,MYJE2 DO 100 I=MYIS,MYIE IF (HBM2(I,J) < 0.5) GO TO 100 ! Ignore columns near lateral boundaries LSFC=LMH(I,J) ! "L" of surface PDSL=PD(I,J)*RES(I,J) ! (Psfc-Ptop)/ETA_sfc CONST=PDSL/GRAV ! (Psfc-Ptop)/(g*ETA_sfc), used for THICK below !--- Initialize column data (1D arrays) IF (CWM(I,J,1) <= CLIMIT) CWM(I,J,1)=EPSQ F_ice(1,I,J)=1. F_rain(1,I,J)=0. F_RimeF(1,I,J)=1. ! NICIALIZACAO DA VARIAVEL T_col T_col(1:LM)=0. DO L=1,LSFC !--- Pressure (Pa) = (Psfc-Ptop)*(ETA/ETA_sfc)+Ptop P_col(L)=PDSL*AETA(L)+PT !--- Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc) THICK_col(L)=CONST*DETA(L) T_col(L)=T(I,J,L) TC=T_col(L)-T0C QV_col(L)=max(EPSQ, Q(I,J,L)) IF (CWM(I,J,L) <= CLIMIT) THEN WC_col(L)=0. IF (TC < T_ICE) THEN F_ice(L,I,J)=1. ELSE F_ice(L,I,J)=0. ENDIF F_rain(L,I,J)=0. F_RimeF(L,I,J)=1. ELSE WC_col(L)=CWM(I,J,L) ENDIF !--- Determine composition of condensate in terms of ! cloud water, ice, & rain WC=WC_col(L) QI=0. QR=0. QW=0. Fice=F_ice(L,I,J) Frain=F_rain(L,I,J) !--- REAL*4 array storage IF (Fice >= 1.) THEN QI=WC ELSE IF (Fice <= 0.) THEN QW=WC ELSE QI=Fice*WC QW=WC-QI ENDIF IF (QW > 0. .AND. Frain > 0.) THEN IF (Frain >= 1.) THEN QR=QW QW=0. ELSE QR=Frain*QW QW=QW-QR ENDIF ENDIF RimeF_col(L)=F_RimeF(L,I,J) ! (real) !! ! !! !--- INTEGER*2 array storage !! ! !! IF (TC.LE.T_ICE .OR. Fice.GE.IFMAX) THEN !! QI=WC !! ELSE IF (Fice .LE. IFMIN) THEN !! QW=WC !! ELSE !! QI=RFSCALE*Fice*WC !! QW=WC-QI !! ENDIF !! IF (QW.GT.0. .AND. Frain.GT.0) THEN !! IF (Frain .EQ. IFMAX) THEN !! QR=QW !! QW=0. !! ELSE !! QR=RFSCALE*Frain*QW !! QW=QW-QR !! ENDIF !! ENDIF !! RimeF_col(L)=F_RimeF(L,I,J)/RIMESCALE QI_col(L)=QI QR_col(L)=QR QW_col(L)=QW ENDDO !####################################################################### !--- Perform the microphysical calculations in this column I_index=I J_index=J CALL GSMCOLUMN ( ARAIN, ASNOW, DTPH, I_index, J_index, LSFC, & P_col, QI_col, QR_col, QV_col, QW_col, RimeF_col, T_col, & THICK_col, WC_col ) !####################################################################### !--- Update storage arrays DO L=1,LSFC TRAIN(I,J,L)=(T_col(L)-T(I,J,L))/DTPH T(I,J,L)=T_col(L) Q(I,J,L)=QV_col(L) TLATGS(I,J,L)=T_col(L)-T(I,J,L) CWM(I,J,L)=WC_col(L) !--- REAL*4 array storage F_RimeF(L,I,J)=MAX(1., RimeF_col(L)) IF (QI_col(L) <= CLIMIT) THEN F_ice(L,I,J)=0. IF (T_col(L) < T_ICEK) F_ice(L,I,J)=1. ELSE F_ice(L,I,J)=MAX( 0., MIN(1., QI_col(L)/WC_col(L)) ) ENDIF IF (QR_col(L) <= CLIMIT) THEN DUM=0 ELSE DUM=QR_col(L)/(QR_col(L)+QW_col(L)) ENDIF F_rain(L,I,J)=DUM !! ! !! !--- INTEGER*2 array storage !! ! !! IF (RimeF_col(L) .LE. 1.) THEN !! I2DUM=IRimeF !! I2DUM=IRimeF !! ELSE !! I2DUM=(RimeF_col(L)+RIMEROUND)*RIMESCALE !! ENDIF !! F_RimeF(L,I,J)=I2DUM !! IF (QI_col(L) .LE. CLIMIT) THEN !! I2DUM=0 !! ELSE !! I2DUM=FSCALE*QI_col(L)/WC_col(L) !! I2DUM=MAX( IFMIN, MIN(IFMAX, I2DUM) ) !! ENDIF !! F_ice(L,I,J)=I2DUM !! IF (QR_col(L) .LE. CLIMIT) THEN !! I2DUM=0 !! ELSE !! I2DUM=FSCALE*QR_col(L)/(QR_col(L)+QW_col(L)) !! I2DUM=MAX( IFMIN, MIN(IFMAX, I2DUM) ) !! ENDIF !! F_rain(L,I,J)=I2DUM ENDDO !--- Update accumulated precipitation statistics !--- Surface precipitation statistics; SR is fraction of surface ! precipitation (if >0) associated with snow APREC(I,J)=(ARAIN+ASNOW)*RRHOL ! Accumulated surface precip (depth in m) !<--- Ying PREC(I,J)=PREC(I,J)+APREC(I,J) ACPREC(I,J)=ACPREC(I,J)+APREC(I,J) IF(APREC(I,J) < 1.E-8) THEN SR(I,J)=0. ELSE SR(I,J)=RRHOL*ASNOW/APREC(I,J) ENDIF !--- Debug statistics IF (PRINT_diag) THEN PRECtot(1)=PRECtot(1)+ARAIN PRECtot(2)=PRECtot(2)+ASNOW PRECmax(1)=MAX(PRECmax(1), ARAIN) PRECmax(2)=MAX(PRECmax(2), ASNOW) ENDIF !####################################################################### !####################################################################### CONTINUE ! End "I" & "J" loops 100 END DO !####################################################################### !####################################################################### !----------------------------------------------------------------------- !--------------------- END of main microphysics loop ------------------- !----------------------------------------------------------------------- time_model=float(NTSD-1)*DT/3600. IF (PRINT_diag .AND. time_model >= Thour_print) THEN CALL MPI_REDUCE(NSTATS,NSTATS_0,ITHILO_N,MPI_INTEGER,MPI_SUM,0, & MPI_COMM_COMP,IRTN) CALL MPI_REDUCE(QMAX,QMAX_0,ITHILO_QM,MPI_REAL,MPI_MAX,0, & MPI_COMM_COMP,IRTN) CALL MPI_REDUCE(PRECmax,PRECmax_0,2,MPI_REAL,MPI_MAX,0, & MPI_COMM_COMP,IRTN) CALL MPI_REDUCE(QTOT,QTOT_0,ITHILO_QT,MPI_REAL,MPI_SUM,0, & MPI_COMM_COMP,IRTN) CALL MPI_REDUCE(PRECtot,PRECtot_0,2,MPI_REAL,MPI_SUM,0, & MPI_COMM_COMP,IRTN) IF (MYPE == 0) THEN HDTPH=3600./DTPH ! Convert precip rates to mm/h DO K=ITLO,ITHI QMAX_0(K,1)=1000.*QMAX_0(K,1) QMAX_0(K,2)=1000.*QMAX_0(K,2) QMAX_0(K,3)=1000.*QMAX_0(K,3) QMAX_0(K,4)=HDTPH*QMAX_0(K,4) QMAX_0(K,5)=HDTPH*QMAX_0(K,5) ENDDO PRECmax_0(1)=HDTPH*PRECmax_0(1) PRECmax_0(2)=HDTPH*PRECmax_0(2) WRITE(6,"(A,F5.2,4(A,G11.4))") '{ Time(h)=',time_model, & ' TRAIN_sfc=',PRECtot_0(1),' TSNOW_sfc=',PRECtot_0(2), & ' RRmax_sfc(mm/h)=',PRECmax_0(1), & ' SRmax_sfc(mm/h)=',PRECmax_0(2) WRITE(6,"(3A)") '{ (C) <--------- Counts ----------> ', & '<----------- g/kg ----------> <----- mm/h ------>', & ' <---- kg/m**2 * # grids ---->' WRITE(6,"(3A)") '{ T NCICE NCMIX NCWAT NCRAIN ', & 'QIMAX QWMAX QRMAX SRMAX RRMAX QITOT ', & 'QWTOT QRTOT' DO K=ITLO,ITHI WRITE(6,"(A,I3,I9,3I7,8G10.4)") & '{ ',K,(NSTATS_0(K,II), II=1,4), & (QMAX_0(K,JJ), JJ=1,5),(QTOT_0(K,KK), KK=1,3) ENDDO WRITE(6,"(3A)") & '{ T TCOND TICND TIEVP TIDEP TREVP ', & 'TRAUT TRACW TIMLT TIACW TIACWI TIACWR ', & 'TIACR' DO K=ITLO,ITHI WRITE(6,"(A,I3,12G10.4)") '{ ',K,(QTOT_0(K,II), II=4,15) ENDDO WRITE(6,"(2A)") & '{ T DEL_QT TVDIF DEL_HYD TWDIF TIDIF ', & 'TRDIF DARAIN DASNOW RimeF' DO K=ITLO,ITHI DEL_HYD=0. DO II=17,19 DEL_HYD=DEL_HYD+QTOT_0(K,II) ENDDO DEL_QT=0. DO II=16,21 DEL_QT=DEL_QT+QTOT_0(K,II) ENDDO IF (QTOT_0(K,22) > 0.) THEN RimeF_bulk=QTOT_0(K,1)/QTOT_0(K,22) ELSE RimeF_bulk=1. ENDIF WRITE(6,"(A,I3,9G10.4)") '{ ',K,DEL_QT,QTOT_0(K,16), & DEL_HYD,(QTOT_0(K,II), II=17,21),RimeF_bulk ENDDO ENDIF NSTATS=0 ! Microphysical statistics dealing w/ grid-point counts QMAX=0. ! Microphysical statistics dealing w/ hydrometeor mass QTOT=0. ! Microphysical statistics dealing w/ hydrometeor mass PRECmax=0. ! Maximum precip rates (rain, snow) at surface (mm/h) PRECtot=0. ! Total precipitation (rain, snow) accumulation at surface Thour_print=Thour_print+DThour_print ENDIF !----------------------------------------------------------------------- !------------------------ Return to main program ----------------------- !----------------------------------------------------------------------- RETURN !----------------------------------------------------------------------- 200 format(a2,i5,f6.2,4(1x,a10,g11.4)) 210 format(a2,i5,f6.2,4(1x,a10,i7)) !----------------------------------------------------------------------- END SUBROUTINE GSMDRIVE