Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.F90 (revision 4665) +++ LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.F90 (revision 4665) @@ -0,0 +1,456 @@ +MODULE lmdz_lscp_tools + + 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 lmdz_lscp_ini, only: iflag_vice, ffallv_con, ffallv_lsc + use lmdz_lscp_ini, 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, temp_cltop, 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 lmdz_lscp_ini, ONLY: t_glace_min, t_glace_max, exposant_glace, iflag_t_glace + USE lmdz_lscp_ini, ONLY : RTT, dist_liq, temp_nowater + + 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 + REAL, INTENT(IN), DIMENSION(klon) :: temp_cltop ! temperature of 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 + + ! if temperature of cloud top <-40°C, + IF (iflag_t_glace .GE. 4) THEN + IF ((temp_cltop(i) .LE. temp_nowater) .AND. (temp(i) .LE. t_glace_max)) THEN + icefrac(i) = 1. + dicefracdT(i) = 0. + 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 lmdz_lscp_ini, 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,temp_cltop) +!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + USE lmdz_lscp_ini, 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, INTENT(OUT), DIMENSION(klon) :: temp_cltop ! temperature of 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 + + DO i=1,klon + k_top = k + 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 + temp_cltop(i) = temp(i,k_top) + ENDDO ! klon + +END SUBROUTINE DISTANCE_TO_CLOUD_TOP +!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +END MODULE lmdz_lscp_tools + + Index: LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90 (revision 4664) +++ (revision ) @@ -1,456 +1,0 @@ -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, temp_cltop, 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, temp_nowater - - 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 - REAL, INTENT(IN), DIMENSION(klon) :: temp_cltop ! temperature of 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 - - ! if temperature of cloud top <-40°C, - IF (iflag_t_glace .GE. 4) THEN - IF ((temp_cltop(i) .LE. temp_nowater) .AND. (temp(i) .LE. t_glace_max)) THEN - icefrac(i) = 1. - dicefracdT(i) = 0. - 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,temp_cltop) -!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - 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, INTENT(OUT), DIMENSION(klon) :: temp_cltop ! temperature of 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 - - DO i=1,klon - k_top = k - 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 - temp_cltop(i) = temp(i,k_top) - ENDDO ! klon - -END SUBROUTINE DISTANCE_TO_CLOUD_TOP -!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - -END MODULE LSCP_TOOLS_MOD - -