MODULE LSCP_TOOLS_MOD IMPLICIT NONE CONTAINS !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE FALLICE_VELOCITY(klon,iwc,temp,rho,pres,ptconv,velo) ! Ref: ! Stubenrauch, C. J., Bonazzola, M., ! Protopapadaki, S. E., & Musat, I. (2019). ! New cloud system metrics to assess bulk ! ice cloud schemes in a GCM. Journal of ! Advances in Modeling Earth Systems, 11, ! 3212–3234. https://doi.org/10.1029/2019MS001642 use lscp_ini_mod, only: iflag_vice, ffallv_con, ffallv_lsc use lscp_ini_mod, only: cice_velo, dice_velo IMPLICIT NONE INTEGER, INTENT(IN) :: klon REAL, INTENT(IN), DIMENSION(klon) :: iwc ! specific ice water content [kg/m3] REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature [K] REAL, INTENT(IN), DIMENSION(klon) :: rho ! dry air density [kg/m3] REAL, INTENT(IN), DIMENSION(klon) :: pres ! air pressure [Pa] LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv ! convective point [-] REAL, INTENT(OUT), DIMENSION(klon) :: velo ! fallspeed velocity of crystals [m/s] INTEGER i REAL logvm,iwcg,tempc,phpa,fallv_tun REAL m2ice, m2snow, vmice, vmsnow REAL aice, bice, asnow, bsnow DO i=1,klon IF (ptconv(i)) THEN fallv_tun=ffallv_con ELSE fallv_tun=ffallv_lsc ENDIF tempc=temp(i)-273.15 ! celcius temp iwcg=MAX(iwc(i)*1000.,1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0 phpa=pres(i)/100. ! pressure in hPa IF (iflag_vice .EQ. 1) THEN ! so-called 'empirical parameterization' in Stubenrauch et al. 2019 if (tempc .GE. -60.0) then logvm= -0.0000414122*tempc*tempc*log(iwcg)-0.00538922*tempc*log(iwcg) & -0.0516344*log(iwcg)+0.00216078*tempc + 1.9714 velo(i)=exp(logvm) else velo(i)=65.0*(iwcg**0.2)*(150./phpa)**0.15 endif velo(i)=fallv_tun*velo(i)/100.0 ! from cm/s to m/s ELSE IF (iflag_vice .EQ. 2) THEN ! so called PSDM empirical coherent bulk ice scheme in Stubenrauch et al. 2019 aice=0.587 bice=2.45 asnow=0.0444 bsnow=2.1 m2ice=((iwcg*0.001/aice)/(exp(13.6-bice*7.76+0.479*bice**2)* & exp((-0.0361+bice*0.0151+0.00149*bice**2)*tempc))) & **(1./(0.807+bice*0.00581+0.0457*bice**2)) vmice=100.*1042.4*exp(13.6-(bice+1)*7.76+0.479*(bice+1.)**2)*exp((-0.0361+ & (bice+1.)*0.0151+0.00149*(bice+1.)**2)*tempc) & *(m2ice**(0.807+(bice+1.)*0.00581+0.0457*(bice+1.)**2))/(iwcg*0.001/aice) vmice=vmice*((1000./phpa)**0.2) m2snow=((iwcg*0.001/asnow)/(exp(13.6-bsnow*7.76+0.479*bsnow**2)* & exp((-0.0361+bsnow*0.0151+0.00149*bsnow**2)*tempc))) & **(1./(0.807+bsnow*0.00581+0.0457*bsnow**2)) vmsnow=100.*14.3*exp(13.6-(bsnow+.416)*7.76+0.479*(bsnow+.416)**2)& *exp((-0.0361+(bsnow+.416)*0.0151+0.00149*(bsnow+.416)**2)*tempc)& *(m2snow**(0.807+(bsnow+.416)*0.00581+0.0457*(bsnow+.416)**2))/(iwcg*0.001/asnow) vmsnow=vmsnow*((1000./phpa)**0.35) velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s ELSE ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990 velo(i) = fallv_tun*cice_velo*((iwcg/1000.)**dice_velo) ENDIF ENDDO END SUBROUTINE FALLICE_VELOCITY !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, icefrac, dicefracdT) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Compute the ice fraction 1-xliq (see e.g. ! Doutriaux-Boucher & Quaas 2004, section 2.2.) ! as a function of temperature ! see also Fig 3 of Madeleine et al. 2020, JAMES !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ USE print_control_mod, ONLY: lunout, prt_level USE lscp_ini_mod, ONLY: t_glace_min, t_glace_max, exposant_glace, iflag_t_glace USE lscp_ini_mod, ONLY : RTT, dist_liq IMPLICIT NONE INTEGER, INTENT(IN) :: klon ! number of horizontal grid points REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature REAL, INTENT(IN), DIMENSION(klon) :: distcltop ! distance to cloud top INTEGER, INTENT(IN) :: iflag_ice_thermo REAL, INTENT(OUT), DIMENSION(klon) :: icefrac REAL, INTENT(OUT), DIMENSION(klon) :: dicefracdT INTEGER i REAL liqfrac_tmp, dicefrac_tmp REAL Dv, denomdep,beta,qsi,dqsidt LOGICAL ice_thermo CHARACTER (len = 20) :: modname = 'lscp_tools' CHARACTER (len = 80) :: abort_message IF ((iflag_t_glace.LT.2) .OR. (iflag_t_glace.GT.6)) THEN abort_message = 'lscp cannot be used if iflag_t_glace<2 or >6' CALL abort_physic(modname,abort_message,1) ENDIF IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN abort_message = 'lscp cannot be used without ice thermodynamics' CALL abort_physic(modname,abort_message,1) ENDIF DO i=1,klon ! old function with sole dependence upon temperature IF (iflag_t_glace .EQ. 2) THEN liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) icefrac(i) = (1.0-liqfrac_tmp)**exposant_glace IF (icefrac(i) .GT.0.) THEN dicefracdT(i)= exposant_glace * (icefrac(i)**(exposant_glace-1.)) & / (t_glace_min - t_glace_max) ENDIF IF ((icefrac(i).EQ.0).OR.(icefrac(i).EQ.1)) THEN dicefracdT(i)=0. ENDIF ENDIF ! function of temperature used in CMIP6 physics IF (iflag_t_glace .EQ. 3) THEN liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) icefrac(i) = 1.0-liqfrac_tmp**exposant_glace IF ((icefrac(i) .GT.0.) .AND. (liqfrac_tmp .GT. 0.)) THEN dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) & / (t_glace_min - t_glace_max) ELSE dicefracdT(i)=0. ENDIF ENDIF ! for iflag_t_glace .GE. 4, the liquid fraction depends upon temperature at cloud top ! and then decreases with decreasing height !with linear function of temperature at cloud top IF (iflag_t_glace .EQ. 4) THEN liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.) dicefrac_tmp = - temp(i)/(t_glace_max-t_glace_min) dicefracdT(i) = dicefrac_tmp*exp(-distcltop(i)/dist_liq) IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN dicefracdT(i) = 0. ENDIF ENDIF ! with CMIP6 function of temperature at cloud top IF (iflag_t_glace .EQ. 5) THEN liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) liqfrac_tmp = liqfrac_tmp**exposant_glace icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.) IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN dicefracdT(i) = 0. ELSE dicefracdT(i) = exposant_glace*((liqfrac_tmp)**(exposant_glace-1.))/(t_glace_min- t_glace_max) & *exp(-distcltop(i)/dist_liq) ENDIF ENDIF ! with modified function of temperature at cloud top ! to get largere values around 260 K, works well with t_glace_min = 241K IF (iflag_t_glace .EQ. 6) THEN IF (temp(i) .GT. t_glace_max) THEN liqfrac_tmp = 1. ELSE liqfrac_tmp = -((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))**2+1. ENDIF liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.) IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN dicefracdT(i) = 0. ELSE dicefracdT(i) = 2*((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))/(t_glace_max-t_glace_min) & *exp(-distcltop(i)/dist_liq) ENDIF ENDIF ENDDO ! klon RETURN END SUBROUTINE ICEFRAC_LSCP !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Calculate qsat following ECMWF method !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ IMPLICIT NONE include "YOMCST.h" include "YOETHF.h" include "FCTTRE.h" INTEGER, INTENT(IN) :: klon ! number of horizontal grid points REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa REAL, INTENT(IN) :: tref ! reference temperature in K LOGICAL, INTENT(IN) :: flagth ! flag for qsat calculation for thermals INTEGER, INTENT(IN) :: phase ! phase: 0=depend on temperature sign (temp>tref -> liquid, tempgammasat*qsat ! Etienne Vignon, March 2021 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ use lscp_ini_mod, only: iflag_gammasat, t_glace_min, RTT IMPLICIT NONE INTEGER, INTENT(IN) :: klon ! number of horizontal grid points REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa REAL, INTENT(OUT), DIMENSION(klon) :: gammasat ! coefficient to multiply qsat with to calculate saturation REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature REAL, DIMENSION(klon) :: qsi,qsl,dqsl,dqsi REAL fcirrus, fac REAL, PARAMETER :: acirrus=2.349 REAL, PARAMETER :: bcirrus=259.0 INTEGER i CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl) CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi) DO i=1,klon IF (temp(i) .GE. RTT) THEN ! warm clouds: condensation at saturation wrt liquid gammasat(i)=1. dgammasatdt(i)=0. ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. t_glace_min)) THEN IF (iflag_gammasat .GE. 2) THEN gammasat(i)=qsl(i)/qsi(i) dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i) ELSE gammasat(i)=1. dgammasatdt(i)=0. ENDIF ELSE IF (iflag_gammasat .GE.1) THEN ! homogeneous freezing of aerosols, according to ! Koop, 2000 and Karcher 2008, QJRMS ! 'Cirrus regime' fcirrus=acirrus-temp(i)/bcirrus IF (fcirrus .LT. qsl(i)/qsi(i)) THEN gammasat(i)=qsl(i)/qsi(i) dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i) ELSE gammasat(i)=fcirrus dgammasatdt(i)=-1.0/bcirrus ENDIF ELSE gammasat(i)=1. dgammasatdt(i)=0. ENDIF ENDIF END DO END SUBROUTINE CALC_GAMMASAT !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D) !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ USE lscp_ini_mod, ONLY : rd,rg,tresh_cl IMPLICIT NONE INTEGER, INTENT(IN) :: klon,klev !number of horizontal and vertical grid points INTEGER, INTENT(IN) :: k ! vertical index REAL, INTENT(IN), DIMENSION(klon,klev) :: temp ! temperature in K REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! pressure middle layer in Pa REAL, INTENT(IN), DIMENSION(klon,klev+1) :: paprs ! pressure interfaces in Pa REAL, INTENT(IN), DIMENSION(klon,klev) :: rneb ! cloud fraction REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D ! distance from cloud top REAL dzlay(klon,klev) REAL zlay(klon,klev) REAL dzinterf INTEGER i,k_top, kvert LOGICAL bool_cl DO i=1,klon ! Initialization height middle of first layer dzlay(i,1) = Rd * temp(i,1) / rg * log(paprs(i,1)/paprs(i,2)) zlay(i,1) = dzlay(i,1)/2 DO kvert=2,klev IF (kvert.EQ.klev) THEN dzlay(i,kvert) = 2*(rd * temp(i,kvert) / rg * log(paprs(i,kvert)/pplay(i,kvert))) ELSE dzlay(i,kvert) = rd * temp(i,kvert) / rg * log(paprs(i,kvert)/paprs(i,kvert+1)) ENDIF dzinterf = rd * temp(i,kvert) / rg * log(pplay(i,kvert-1)/pplay(i,kvert)) zlay(i,kvert) = zlay(i,kvert-1) + dzinterf ENDDO ENDDO k_top = k DO i=1,klon IF (rneb(i,k) .LE. tresh_cl) THEN bool_cl = .FALSE. ELSE bool_cl = .TRUE. ENDIF DO WHILE ((bool_cl) .AND. (k_top .LE. klev)) ! find cloud top IF (rneb(i,k_top) .GT. tresh_cl) THEN k_top = k_top + 1 ELSE bool_cl = .FALSE. k_top = k_top - 1 ENDIF ENDDO k_top=min(k_top,klev) !dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to !interf for layer of cloud top distcltop1D(i) = zlay(i,k_top) - zlay(i,k) + dzlay(i,k_top)/2 ENDDO ! klon END SUBROUTINE DISTANCE_TO_CLOUD_TOP !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ END MODULE LSCP_TOOLS_MOD