Changeset 4559
- Timestamp:
- Jun 6, 2023, 12:09:27 PM (12 months ago)
- Location:
- LMDZ6/trunk/libf/phylmd
- Files:
-
- 3 edited
-
lscp_ini_mod.F90 (modified) (13 diffs)
-
lscp_mod.F90 (modified) (10 diffs)
-
lscp_tools_mod.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/lscp_ini_mod.F90
r4535 r4559 12 12 !$OMP THREADPRIVATE(seuil_neb) 13 13 14 INTEGER, SAVE :: ni nter=5! number of iterations to calculate autoconversion to precipitation15 !$OMP THREADPRIVATE(ni nter)14 INTEGER, SAVE :: niter_lscp=5 ! number of iterations to calculate autoconversion to precipitation 15 !$OMP THREADPRIVATE(niter_lscp) 16 16 17 17 INTEGER,SAVE :: iflag_evap_prec=1 ! precipitation evaporation flag. 0: nothing, 1: "old way", … … 41 41 !$OMP THREADPRIVATE(ok_radocond_snow) 42 42 43 LOGICAL, SAVE :: ok_debug_autoconversion=.true. ! removes a bug in the autoconversion process44 !$OMP THREADPRIVATE(ok_debug_autoconversion)45 46 43 REAL, SAVE :: t_glace_min=258.0 ! lower-bound temperature parameter for cloud phase determination 47 44 !$OMP THREADPRIVATE(t_glace_min) … … 77 74 !$OMP THREADPRIVATE(iflag_pdf) 78 75 76 INTEGER, SAVE :: iflag_autoconversion=0 ! autoconversion option 77 !$OMP THREADPRIVATE(iflag_autoconversion) 78 79 79 LOGICAL, SAVE :: reevap_ice=.false. ! no liquid precip for T< threshold 80 80 !$OMP THREADPRIVATE(reevap_ice) … … 92 92 !$OMP THREADPRIVATE(cld_tau_con) 93 93 94 REAL, SAVE :: cld_expo_lsc=2. ! liquid autoconversion threshold exponent, stratiform rain 95 !$OMP THREADPRIVATE(cld_expo_lsc) 96 97 REAL, SAVE :: cld_expo_con=2. ! liquid autoconversion threshold exponent, convective rain 98 !$OMP THREADPRIVATE(cld_expo_con) 99 94 100 REAL, SAVE :: ffallv_lsc=1. ! tuning coefficient crystal fall velocity, stratiform 95 101 !$OMP THREADPRIVATE(ffallv_lsc) … … 104 110 !$OMP THREADPRIVATE(coef_eva_i) 105 111 106 112 REAL cice_velo ! factor in the ice fall velocity formulation 113 PARAMETER (cice_velo=1.645) 114 115 REAL dice_velo ! exponent in the ice fall velocity formulation 116 PARAMETER (dice_velo=0.16) 107 117 108 118 … … 122 132 REAL, INTENT(IN) :: RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in 123 133 REAL, INTENT(IN) :: RVTMP2_in, RTT_in, RD_in, RG_in 134 character (len=20) :: modname='lscp_ini_mod' 135 character (len=80) :: abort_message 124 136 125 137 … … 136 148 137 149 138 CALL getin_p('ni nter',ninter)150 CALL getin_p('niter_lscp',niter_lscp) 139 151 CALL getin_p('iflag_evap_prec',iflag_evap_prec) 140 152 CALL getin_p('seuil_neb',seuil_neb) … … 142 154 CALL getin_p('iflag_mpc_bl',iflag_mpc_bl) 143 155 CALL getin_p('ok_radocond_snow',ok_radocond_snow) 144 CALL getin_p('ok_debug_autoconversion',ok_debug_autoconversion)145 156 CALL getin_p('t_glace_max',t_glace_max) 146 157 CALL getin_p('t_glace_min',t_glace_min) … … 159 170 CALL getin_p('cld_tau_lsc',cld_tau_lsc) 160 171 CALL getin_p('cld_tau_con',cld_tau_con) 172 CALL getin_p('cld_expo_lsc',cld_expo_lsc) 173 CALL getin_p('cld_expo_con',cld_expo_con) 161 174 CALL getin_p('ffallv_lsc',ffallv_lsc) 162 175 CALL getin_p('ffallv_lsc',ffallv_con) … … 164 177 coef_eva_i=coef_eva 165 178 CALL getin_p('coef_eva_i',coef_eva_i) 166 167 168 169 170 WRITE(lunout,*) 'lscp, ninter:', ninter 179 CALL getin_p('iflag_autoconversion',iflag_autoconversion) 180 181 182 183 184 WRITE(lunout,*) 'lscp, niter_lscp:', niter_lscp 171 185 WRITE(lunout,*) 'lscp, iflag_evap_prec:', iflag_evap_prec 172 186 WRITE(lunout,*) 'lscp, seuil_neb:', seuil_neb … … 174 188 WRITE(lunout,*) 'lscp, iflag_mpc_bl:', iflag_mpc_bl 175 189 WRITE(lunout,*) 'lscp, ok_radocond_snow:', ok_radocond_snow 176 WRITE(lunout,*) 'lscp, ok_debug_autoconversion:', ok_debug_autoconversion177 190 WRITE(lunout,*) 'lscp, t_glace_max:', t_glace_max 178 191 WRITE(lunout,*) 'lscp, t_glace_min:', t_glace_min … … 191 204 WRITE(lunout,*) 'lscp, cld_tau_lsc', cld_tau_lsc 192 205 WRITE(lunout,*) 'lscp, cld_tau_con', cld_tau_con 206 WRITE(lunout,*) 'lscp, cld_expo_lsc', cld_expo_lsc 207 WRITE(lunout,*) 'lscp, cld_expo_con', cld_expo_con 193 208 WRITE(lunout,*) 'lscp, ffallv_lsc', ffallv_lsc 194 209 WRITE(lunout,*) 'lscp, ffallv_con', ffallv_con 195 210 WRITE(lunout,*) 'lscp, coef_eva', coef_eva 196 211 WRITE(lunout,*) 'lscp, coef_eva_i', coef_eva_i 212 WRITE(lunout,*) 'lscp, iflag_autoconversion', iflag_autoconversion 197 213 198 214 … … 201 217 202 218 ! check for precipitation sub-time steps 203 IF (ABS(dtime/REAL(ni nter)-360.0).GT.0.001) THEN219 IF (ABS(dtime/REAL(niter_lscp)-360.0).GT.0.001) THEN 204 220 WRITE(lunout,*) 'lscp: it is not expected, see Z.X.Li', dtime 205 221 WRITE(lunout,*) 'I would prefer a 6 min sub-timestep' 222 ENDIF 223 224 ! check consistency between numerical resolution of autoconversion 225 ! and other options 226 227 IF (iflag_autoconversion .EQ. 2) THEN 228 IF ((iflag_vice .NE. 0) .OR. (niter_lscp .GT. 1)) THEN 229 abort_message = 'in lscp, iflag_autoconversion=2 requires iflag_vice=0 and niter_lscp=1' 230 CALL abort_physic (modname,abort_message,1) 231 ENDIF 206 232 ENDIF 207 233 -
LMDZ6/trunk/libf/phylmd/lscp_mod.F90
r4535 r4559 93 93 USE ice_sursat_mod, ONLY : ice_sursat 94 94 95 USE lscp_ini_mod, ONLY : seuil_neb, ni nter, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min96 USE lscp_ini_mod, ONLY : iflag_mpc_bl, ok_radocond_snow, a_tr_sca, ok_debug_autoconversion95 USE lscp_ini_mod, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min 96 USE lscp_ini_mod, ONLY : iflag_mpc_bl, ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc 97 97 USE lscp_ini_mod, ONLY : iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min 98 98 USE lscp_ini_mod, ONLY : coef_eva, coef_eva_i,cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con 99 USE lscp_ini_mod, ONLY : iflag_bergeron, iflag_fisrtilp_qsat 99 USE lscp_ini_mod, ONLY : iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo 100 USE lscp_ini_mod, ONLY : iflag_autoconversion, ffallv_con, ffallv_lsc 100 101 USE lscp_ini_mod, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG 101 102 … … 190 191 191 192 REAL qsl(klon), qsi(klon) 192 REAL zct, zcl 193 REAL zct, zcl,zexpo 193 194 REAL ctot(klon,klev) 194 195 REAL ctot_vol(klon,klev) … … 232 233 REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon) 233 234 REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon) 234 REAL velo(klon,klev), vr 235 REAL velo(klon,klev), vr, ffallv 235 236 REAL qlmpc, qimpc, rnebmpc 236 237 REAL radocondi(klon,klev), radocondl(klon,klev) 238 REAL effective_zneb 237 239 238 240 INTEGER i, k, n, kk … … 922 924 ENDIF 923 925 926 927 ! Autoconversion 928 ! ------------------------------------------------------------------------------- 929 930 924 931 ! Initialisation of zoliq and radocond variables 925 932 … … 933 940 iwc(i) = 0. 934 941 zneb(i) = MAX(rneb(i,k),seuil_neb) 935 radocond(i,k) = zoliq(i)/REAL(ni nter+1)942 radocond(i,k) = zoliq(i)/REAL(niter_lscp+1) 936 943 radicefrac(i,k) = zfice(i) 937 radocondi(i,k)=zoliq(i)*zfice(i)/REAL(ni nter+1)938 radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(ni nter+1)944 radocondi(i,k)=zoliq(i)*zfice(i)/REAL(niter_lscp+1) 945 radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1) 939 946 ENDDO 940 947 941 ! Autoconversion 942 ! ------------------------------------------------------------------------------- 943 944 DO n = 1, ninter 948 949 DO n = 1, niter_lscp 945 950 946 951 DO i=1,klon … … 968 973 ELSE 969 974 970 IF (ok_debug_autoconversion) THEN ! if condition to be removed after test phase971 972 ! water quantity to remove: zchau (Sundqvist, 1978)973 ! same thing for the ice: zfroi (Zender & Kiehl, 1997)974 975 IF (ptconv(i,k)) THEN ! if convective point 975 976 zcl=cld_lc_con 976 977 zct=1./cld_tau_con 977 zfroi = dtime/REAL(ninter)/zdz(i)*zoliqi(i)*velo(i,k) 978 zexpo=cld_expo_con 979 ffallv=ffallv_con 978 980 ELSE 979 981 zcl=cld_lc_lsc 980 982 zct=1./cld_tau_lsc 981 z froi = dtime/REAL(ninter)/zdz(i)*zoliqi(i) & ! dqice/dt=1/rho*d(rho*wice*qice)/dz982 *velo(i,k)983 zexpo=cld_expo_lsc 984 ffallv=ffallv_lsc 983 985 ENDIF 986 984 987 985 988 ! if vertical heterogeneity is taken into account, we use … … 988 991 ! reduces the in-cloud water). 989 992 993 ! Liquid water quantity to remove: zchau (Sundqvist, 1978) 994 ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2) 995 !......................................................... 990 996 IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN 991 zchau = zct *dtime/REAL(ninter) * zoliql(i) &992 *(1.0-EXP(-(zoliql(i)/ctot_vol(i,k)/zcl)**2))993 ELSE994 zchau = zct *dtime/REAL(ninter) * zoliql(i) &995 *(1.0-EXP(-(zoliql(i)/zneb(i)/zcl)**2)) ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)996 ENDIF997 998 ELSE ! with old bug in autoconversion999 1000 ! water quantity to remove: zchau (Sundqvist, 1978)1001 ! same thing for the ice: zfroi (Zender & Kiehl, 1997)1002 IF (ptconv(i,k)) THEN ! if convective point1003 zcl=cld_lc_con1004 zct=1./cld_tau_con1005 zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i)*velo(i,k)*zfice(i)1006 ELSE1007 zcl=cld_lc_lsc1008 zct=1./cld_tau_lsc1009 zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) & ! dqice/dt=1/rho*d(rho*wice*qice)/dz1010 *velo(i,k) * zfice(i)1011 ENDIF1012 997 1013 998 ! if vertical heterogeneity is taken into account, we use … … 1015 1000 ! surface fraction (which is larger and artificially 1016 1001 ! reduces the in-cloud water). 1017 1018 1019 IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN 1020 zchau = zct *dtime/REAL(ninter) * zoliq(i) & 1021 *(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl)**2)) *(1.-zfice(i)) 1002 effective_zneb=ctot_vol(i,k) 1022 1003 ELSE 1023 zchau = zct *dtime/REAL(ninter) * zoliq(i) & 1024 *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl)**2)) *(1.-zfice(i)) ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2) 1004 effective_zneb=zneb(i) 1005 ENDIF 1006 1007 1008 IF (iflag_autoconversion .EQ. 2) THEN 1009 ! two-steps resolution with niter_lscp=1 sufficient 1010 ! we first treat the second term (with exponential) in an explicit way 1011 ! and then treat the first term (-q/tau) in an exact way 1012 zchau=zoliql(i)*(1.-exp(-dtime/REAL(niter_lscp)*zct & 1013 *(1.-exp(-(zoliql(i)/effective_zneb/zcl)**zexpo)))) 1014 ELSE 1015 ! old explicit resolution with subtimesteps 1016 zchau = zct*dtime/REAL(niter_lscp)*zoliql(i) & 1017 *(1.0-EXP(-(zoliql(i)/effective_zneb/zcl)**zexpo)) 1025 1018 ENDIF 1026 1019 1027 ENDIF ! ok_debug_autoconversion 1020 1021 ! Ice water quantity to remove (Zender & Kiehl, 1997) 1022 ! dqice/dt=1/rho*d(rho*wice*qice)/dz 1023 !.................................... 1024 IF (iflag_autoconversion .EQ. 2) THEN 1025 ! exact resolution, niter_lscp=1 is sufficient but works only 1026 ! with iflag_vice=0 1027 IF (zoliqi(i) .GT. 0.) THEN 1028 zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) & 1029 +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)*ffallv)**(-1./dice_velo)) 1030 ELSE 1031 zfroi=0. 1032 ENDIF 1033 ELSE 1034 ! old explicit resolution with subtimesteps 1035 zfroi = dtime/REAL(niter_lscp)/zdz(i)*zoliqi(i)*velo(i,k) 1036 ENDIF 1028 1037 1029 1038 zrain = MIN(MAX(zchau,0.0),zoliql(i)) … … 1033 1042 ENDIF 1034 1043 1035 zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain , 0.0) 1036 zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow , 0.0) 1037 zoliq(i) = MAX(zoliq(i)-zprecip , 0.0) 1044 IF (iflag_autoconversion .GE. 1) THEN 1045 ! debugged version with phase conservation through the autoconversion process 1046 zoliql(i) = MAX(zoliql(i)-1.*zrain , 0.0) 1047 zoliqi(i) = MAX(zoliqi(i)-1.*zsnow , 0.0) 1048 zoliq(i) = MAX(zoliq(i)-zprecip , 0.0) 1049 ELSE 1050 ! bugged version with phase resetting 1051 zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain , 0.0) 1052 zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow , 0.0) 1053 zoliq(i) = MAX(zoliq(i)-zprecip , 0.0) 1054 ENDIF 1038 1055 1039 1056 ! c_iso: call isotope_conversion (for readibility) or duplicate 1040 1057 1041 radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(ni nter+1)1042 radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(ni nter+1)1043 radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(ni nter+1)1058 radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(niter_lscp+1) 1059 radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(niter_lscp+1) 1060 radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(niter_lscp+1) 1044 1061 1045 1062 ENDIF ! rneb >0 … … 1137 1154 ENDDO 1138 1155 1156 1139 1157 ENDIF 1140 1158 -
LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90
r4535 r4559 17 17 18 18 use lscp_ini_mod, only: iflag_vice, ffallv_con, ffallv_lsc 19 use lscp_ini_mod, only: cice_velo, dice_velo 19 20 20 21 IMPLICIT NONE … … 31 32 32 33 INTEGER i 33 REAL logvm,iwcg,tempc,phpa, cvel,dvel,fallv_tun34 REAL logvm,iwcg,tempc,phpa,fallv_tun 34 35 REAL m2ice, m2snow, vmice, vmsnow 35 36 REAL aice, bice, asnow, bsnow … … 60 61 61 62 velo(i)=fallv_tun*velo(i)/100.0 ! from cm/s to m/s 62 dvel=0.263 cvel=fallv_tun*65.0*(rho(i)**0.2)*(150./phpa)**0.1564 63 65 64 ELSE IF (iflag_vice .EQ. 2) THEN … … 92 91 vmsnow=vmsnow*((1000./phpa)**0.35) 93 92 velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s 94 dvel=0.295 cvel=velo(i)/((iwcg/1000.*rho(i))**dvel)96 93 97 94 ELSE 98 95 ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990 99 velo(i) = fallv_tun*3.29/2.0*((iwcg/1000.)**0.16) 100 dvel=0.16 101 cvel=fallv_tun*3.29/2.0*(rho(i)**0.16) 96 velo(i) = fallv_tun*cice_velo*((iwcg/1000.)**dice_velo) 102 97 ENDIF 103 98 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
