Changeset 5144 for LMDZ6/branches/Amaury_dev/libf/phylmd/nuage.F90
- Timestamp:
- Jul 29, 2024, 11:01:04 PM (3 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/phylmd/nuage.F90
r5143 r5144 1 1 ! $Id$ 2 2 3 SUBROUTINE nuage(paprs, pplay, t, pqlwp, picefra, pclc, pcltau, pclemi, pch, pcl, pcm, &4 pct, pctlwp, ok_aie, mass_solu_aero, mass_solu_aero_pi, bl95_b0, bl95_b1, distcltop, &5 temp_cltop, cldtaupi, re, fl)3 SUBROUTINE nuage(paprs, pplay, t, pqlwp, picefra, pclc, pcltau, pclemi, pch, pcl, pcm, & 4 pct, pctlwp, ok_aie, mass_solu_aero, mass_solu_aero_pi, bl95_b0, bl95_b1, distcltop, & 5 temp_cltop, cldtaupi, re, fl) 6 6 USE dimphy 7 7 USE lmdz_lscp_tools, ONLY: icefrac_lscp … … 11 11 USE lmdz_clesphys 12 12 USE lmdz_nuage_params ! JBM 3/14 13 USE lmdz_yomcst 13 14 14 15 IMPLICIT NONE … … 41 42 ! ====================================================================== 42 43 43 include "YOMCST.h" 44 45 REAL paprs(klon, klev+1), pplay(klon, klev) 44 REAL paprs(klon, klev + 1), pplay(klon, klev) 46 45 REAL t(klon, klev) 47 46 48 47 REAL pclc(klon, klev) 49 REAL pqlwp(klon, klev), picefra(klon, klev)48 REAL pqlwp(klon, klev), picefra(klon, klev) 50 49 REAL pcltau(klon, klev), pclemi(klon, klev) 51 50 52 51 REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon) 53 REAL distcltop(klon, klev)54 REAL temp_cltop(klon, klev)52 REAL distcltop(klon, klev) 53 REAL temp_cltop(klon, klev) 55 54 LOGICAL lo 56 55 57 56 REAL cetahb, cetamb 58 PARAMETER (cetahb =0.45, cetamb=0.80)57 PARAMETER (cetahb = 0.45, cetamb = 0.80) 59 58 60 59 INTEGER i, k … … 62 61 63 62 REAL radius, rad_chaud 64 ! JBM (3/14) parameters already defined in nuage.h:65 ! REAL rad_froid, rad_chau1, rad_chau266 ! PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)63 ! JBM (3/14) parameters already defined in nuage.h: 64 ! REAL rad_froid, rad_chau1, rad_chau2 65 ! PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0) 67 66 ! cc PARAMETER (rad_chaud=15.0, rad_froid=35.0) 68 67 ! sintex initial PARAMETER (rad_chaud=10.0, rad_froid=30.0) 69 68 REAL coef, coef_froi, coef_chau 70 PARAMETER (coef_chau =0.13, coef_froi=0.09)69 PARAMETER (coef_chau = 0.13, coef_froi = 0.09) 71 70 REAL seuil_neb 72 PARAMETER (seuil_neb =0.001)73 ! JBM (3/14) nexpo is replaced by exposant_glace74 ! REAL nexpo ! exponentiel pour glace/eau75 ! PARAMETER (nexpo=6.)71 PARAMETER (seuil_neb = 0.001) 72 ! JBM (3/14) nexpo is replaced by exposant_glace 73 ! REAL nexpo ! exponentiel pour glace/eau 74 ! PARAMETER (nexpo=6.) 76 75 REAL, PARAMETER :: t_glace_min_old = 258. 77 76 INTEGER, PARAMETER :: exposant_glace_old = 6 … … 104 103 105 104 DO k = 1, klev 106 107 108 zfice(i) = 1.0 - (t(i, k)-t_glace_min_old)/(273.13-t_glace_min_old)109 zfice(i) = min(max(zfice(i), 0.0), 1.0)105 IF (iflag_t_glace==0) THEN 106 DO i = 1, klon 107 zfice(i) = 1.0 - (t(i, k) - t_glace_min_old) / (273.13 - t_glace_min_old) 108 zfice(i) = min(max(zfice(i), 0.0), 1.0) 110 109 zfice(i) = zfice(i)**exposant_glace_old 111 112 113 ! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod114 ! zfice(i) = icefrac_lsc(t(i,k), t_glace_min, &115 ! t_glace_max, exposant_glace)116 117 CALL icefrac_lscp(klon,t(:,k),iflag_ice_thermo,distcltop(:,k),temp_cltop(:,k),zfice(:),dzfice(:))118 119 CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:))120 121 122 123 IF (ok_new_lscp .AND. ok_icefra_lscp) THEN124 ! EV: take the ice fraction directly from the lscp code125 ! consistent only for non convective grid points126 ! critical for mixed phase clouds127 DO i =1,klon128 IF (.NOT. ptconv(i,k)) THEN129 zfice(i)=picefra(i,k)130 ENDIF110 ENDDO 111 ELSE ! of IF (iflag_t_glace.EQ.0) 112 ! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod 113 ! zfice(i) = icefrac_lsc(t(i,k), t_glace_min, & 114 ! t_glace_max, exposant_glace) 115 IF (ok_new_lscp) THEN 116 CALL icefrac_lscp(klon, t(:, k), iflag_ice_thermo, distcltop(:, k), temp_cltop(:, k), zfice(:), dzfice(:)) 117 ELSE 118 CALL icefrac_lsc(klon, t(:, k), pplay(:, k) / paprs(:, 1), zfice(:)) 119 120 ENDIF 121 122 IF (ok_new_lscp .AND. ok_icefra_lscp) THEN 123 ! EV: take the ice fraction directly from the lscp code 124 ! consistent only for non convective grid points 125 ! critical for mixed phase clouds 126 DO i = 1, klon 127 IF (.NOT. ptconv(i, k)) THEN 128 zfice(i) = picefra(i, k) 129 ENDIF 131 130 ENDDO 132 ENDIF 133 134 135 ENDIF 131 ENDIF 132 133 ENDIF 136 134 137 135 DO i = 1, klon … … 139 137 IF (k<=3) rad_chaud = rad_chau2 140 138 141 pclc(i, k) = max(pclc(i, k), seuil_neb)142 zflwp = 1000. *pqlwp(i, k)/rg/pclc(i, k)*(paprs(i,k)-paprs(i,k+1))139 pclc(i, k) = max(pclc(i, k), seuil_neb) 140 zflwp = 1000. * pqlwp(i, k) / rg / pclc(i, k) * (paprs(i, k) - paprs(i, k + 1)) 143 141 144 142 IF (ok_aie) THEN 145 146 147 cdnc(i, k) = 10.**(bl95_b0 +bl95_b1*log(max(mass_solu_aero(i,k), &148 1.E-4))/log(10.))*1.E6 !-m-3149 150 151 152 cdnc(i, k) = min(1000.E6, max(20.E6, cdnc(i,k)))153 cdnc_pi(i, k) = 10.**(bl95_b0 +bl95_b1*log(max(mass_solu_aero_pi(i,k), &154 1.E-4))/log(10.))*1.E6 !-m-3155 cdnc_pi(i, k) = min(1000.E6, max(20.E6, cdnc_pi(i,k)))156 157 158 159 160 161 162 163 rad_chaud = 1.1 *((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi*1000. &164 *cdnc(i,k)))**(1./3.)165 166 167 168 rad_chaud = max(rad_chaud *1.E6, 3.)169 170 171 172 173 174 175 176 177 178 179 180 fl(i, k) = pclc(i, k) *(1.-zfice(i))181 re(i, k) = rad_chaud *fl(i, k)182 183 184 185 186 187 188 radius = max(1.1E6 *((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi* &189 1000.*cdnc_pi(i,k)))**(1./3.), 3.)*(1.-zfice(i)) + rad_froid*zfice(i)190 cldtaupi(i, k) = 3.0 /2.0*zflwp/radius143 ! Formula "D" of Boucher and Lohmann, Tellus, 1995 144 145 cdnc(i, k) = 10.**(bl95_b0 + bl95_b1 * log(max(mass_solu_aero(i, k), & 146 1.E-4)) / log(10.)) * 1.E6 !-m-3 147 ! Cloud droplet number concentration (CDNC) is restricted 148 ! to be within [20, 1000 cm^3] 149 150 cdnc(i, k) = min(1000.E6, max(20.E6, cdnc(i, k))) 151 cdnc_pi(i, k) = 10.**(bl95_b0 + bl95_b1 * log(max(mass_solu_aero_pi(i, k), & 152 1.E-4)) / log(10.)) * 1.E6 !-m-3 153 cdnc_pi(i, k) = min(1000.E6, max(20.E6, cdnc_pi(i, k))) 154 155 156 ! air density: pplay(i,k) / (RD * zT(i,k)) 157 ! factor 1.1: derive effective radius from volume-mean radius 158 ! factor 1000 is the water density 159 ! _chaud means that this is the CDR for liquid water clouds 160 161 rad_chaud = 1.1 * ((pqlwp(i, k) * pplay(i, k) / (rd * t(i, k))) / (4. / 3. * rpi * 1000. & 162 * cdnc(i, k)))**(1. / 3.) 163 164 ! Convert to um. CDR shall be at least 3 um. 165 166 rad_chaud = max(rad_chaud * 1.E6, 3.) 167 168 ! For output diagnostics 169 170 ! Cloud droplet effective radius [um] 171 172 ! we multiply here with f * xl (fraction of liquid water 173 ! clouds in the grid cell) to avoid problems in the 174 ! averaging of the output. 175 ! In the output of IOIPSL, derive the real cloud droplet 176 ! effective radius as re/fl 177 178 fl(i, k) = pclc(i, k) * (1. - zfice(i)) 179 re(i, k) = rad_chaud * fl(i, k) 180 181 ! Pre-industrial cloud opt thickness 182 183 ! "radius" is calculated as rad_chaud above (plus the 184 ! ice cloud contribution) but using cdnc_pi instead of 185 ! cdnc. 186 radius = max(1.1E6 * ((pqlwp(i, k) * pplay(i, k) / (rd * t(i, k))) / (4. / 3. * rpi * & 187 1000. * cdnc_pi(i, k)))**(1. / 3.), 3.) * (1. - zfice(i)) + rad_froid * zfice(i) 188 cldtaupi(i, k) = 3.0 / 2.0 * zflwp / radius 191 189 END IF ! ok_aie 192 190 193 radius = rad_chaud *(1.-zfice(i)) + rad_froid*zfice(i)194 coef = coef_chau *(1.-zfice(i)) + coef_froi*zfice(i)195 pcltau(i, k) = 3.0 /2.0*zflwp/radius196 pclemi(i, k) = 1.0 - exp(-coef *zflwp)197 lo = (pclc(i, k)<=seuil_neb)191 radius = rad_chaud * (1. - zfice(i)) + rad_froid * zfice(i) 192 coef = coef_chau * (1. - zfice(i)) + coef_froi * zfice(i) 193 pcltau(i, k) = 3.0 / 2.0 * zflwp / radius 194 pclemi(i, k) = 1.0 - exp(-coef * zflwp) 195 lo = (pclc(i, k)<=seuil_neb) 198 196 IF (lo) pclc(i, k) = 0.0 199 197 IF (lo) pcltau(i, k) = 0.0 … … 241 239 DO k = klev, 1, -1 242 240 DO i = 1, klon 243 pctlwp(i) = pctlwp(i) + pqlwp(i, k) *(paprs(i,k)-paprs(i,k+1))/rg244 pct(i) = pct(i) *(1.0-pclc(i,k))245 IF (pplay(i, k)<=cetahb*paprs(i,1)) pch(i) = pch(i)*(1.0-pclc(i,k))246 IF (pplay(i, k)>cetahb*paprs(i,1) .AND. pplay(i,k)<=cetamb*paprs(i,1)) &247 pcm(i) = pcm(i)*(1.0-pclc(i,k))248 IF (pplay(i, k)>cetamb*paprs(i,1)) pcl(i) = pcl(i)*(1.0-pclc(i,k))241 pctlwp(i) = pctlwp(i) + pqlwp(i, k) * (paprs(i, k) - paprs(i, k + 1)) / rg 242 pct(i) = pct(i) * (1.0 - pclc(i, k)) 243 IF (pplay(i, k)<=cetahb * paprs(i, 1)) pch(i) = pch(i) * (1.0 - pclc(i, k)) 244 IF (pplay(i, k)>cetahb * paprs(i, 1) .AND. pplay(i, k)<=cetamb * paprs(i, 1)) & 245 pcm(i) = pcm(i) * (1.0 - pclc(i, k)) 246 IF (pplay(i, k)>cetamb * paprs(i, 1)) pcl(i) = pcl(i) * (1.0 - pclc(i, k)) 249 247 END DO 250 248 END DO … … 257 255 END DO 258 256 259 260 257 END SUBROUTINE nuage 261 258 SUBROUTINE diagcld1(paprs, pplay, rain, snow, kbot, ktop, diafra, dialiq) 262 259 USE dimphy 260 USE lmdz_yomcst 261 263 262 IMPLICIT NONE 264 263 … … 270 269 ! ces nuages. Je dois avouer que c'est une frustration. 271 270 272 include "YOMCST.h"273 274 271 ! Arguments d'entree: 275 REAL paprs(klon, klev +1) ! pression (Pa) a inter-couche272 REAL paprs(klon, klev + 1) ! pression (Pa) a inter-couche 276 273 REAL pplay(klon, klev) ! pression (Pa) au milieu de couche 277 274 REAL t(klon, klev) ! temperature (K) … … 288 285 ! Constantes ajustables: 289 286 REAL canva, canvb, canvh 290 PARAMETER (canva =2.0, canvb=0.3, canvh=0.4)287 PARAMETER (canva = 2.0, canvb = 0.3, canvh = 0.4) 291 288 REAL cca, ccb, ccc 292 PARAMETER (cca =0.125, ccb=1.5, ccc=0.8)289 PARAMETER (cca = 0.125, ccb = 1.5, ccc = 0.8) 293 290 REAL ccfct, ccscal 294 PARAMETER (ccfct =0.400)295 PARAMETER (ccscal =1.0E+11)291 PARAMETER (ccfct = 0.400) 292 PARAMETER (ccscal = 1.0E+11) 296 293 REAL cetahb, cetamb 297 PARAMETER (cetahb =0.45, cetamb=0.80)294 PARAMETER (cetahb = 0.45, cetamb = 0.80) 298 295 REAL cclwmr 299 PARAMETER (cclwmr =1.E-04)296 PARAMETER (cclwmr = 1.E-04) 300 297 REAL zepscr 301 PARAMETER (zepscr =1.0E-10)298 PARAMETER (zepscr = 1.0E-10) 302 299 303 300 ! Variables locales: … … 316 313 DO i = 1, klon ! Calculer la fraction nuageuse 317 314 zcc(i) = 0.0 318 IF ((rain(i) +snow(i))>0.) THEN319 zcc(i) = cca *log(max(zepscr,(rain(i)+snow(i))*ccscal)) - ccb320 zcc(i) = min(ccc, max(0.0, zcc(i)))315 IF ((rain(i) + snow(i))>0.) THEN 316 zcc(i) = cca * log(max(zepscr, (rain(i) + snow(i)) * ccscal)) - ccb 317 zcc(i) = min(ccc, max(0.0, zcc(i))) 321 318 END IF 322 319 END DO 323 320 324 321 DO i = 1, klon ! pour traiter les enclumes 325 diafra(i, ktop(i)) = max(diafra(i, ktop(i)), zcc(i)*ccfct)326 IF ((zcc(i)>=canvh) .AND. (pplay(i, ktop(i))<=cetahb*paprs(i, &327 1))) diafra(i, ktop(i)) = max(diafra(i,ktop(i)), max(zcc(&328 i)*ccfct,canva*(zcc(i)-canvb)))329 dialiq(i, ktop(i)) = cclwmr *diafra(i, ktop(i))322 diafra(i, ktop(i)) = max(diafra(i, ktop(i)), zcc(i) * ccfct) 323 IF ((zcc(i)>=canvh) .AND. (pplay(i, ktop(i))<=cetahb * paprs(i, & 324 1))) diafra(i, ktop(i)) = max(diafra(i, ktop(i)), max(zcc(& 325 i) * ccfct, canva * (zcc(i) - canvb))) 326 dialiq(i, ktop(i)) = cclwmr * diafra(i, ktop(i)) 330 327 END DO 331 328 … … 333 330 DO i = 1, klon 334 331 IF (k<ktop(i) .AND. k>=kbot(i)) THEN 335 diafra(i, k) = max(diafra(i, k), zcc(i)*ccfct)336 dialiq(i, k) = cclwmr *diafra(i, k)332 diafra(i, k) = max(diafra(i, k), zcc(i) * ccfct) 333 dialiq(i, k) = cclwmr * diafra(i, k) 337 334 END IF 338 335 END DO 339 336 END DO 340 341 337 342 338 END SUBROUTINE diagcld1 343 339 SUBROUTINE diagcld2(paprs, pplay, t, q, diafra, dialiq) 344 340 USE dimphy 345 USE lmdz_ YOETHF341 USE lmdz_yoethf 346 342 USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep 343 USE lmdz_yomcst 347 344 348 345 IMPLICIT NONE 349 346 350 include "YOMCST.h"351 352 347 ! Arguments d'entree: 353 REAL paprs(klon, klev +1) ! pression (Pa) a inter-couche348 REAL paprs(klon, klev + 1) ! pression (Pa) a inter-couche 354 349 REAL pplay(klon, klev) ! pression (Pa) au milieu de couche 355 350 REAL t(klon, klev) ! temperature (K) … … 361 356 362 357 REAL cetamb 363 PARAMETER (cetamb =0.80)358 PARAMETER (cetamb = 0.80) 364 359 REAL cloia, cloib, cloic, cloid 365 PARAMETER (cloia =1.0E+02, cloib=-10.00, cloic=-0.6, cloid=5.0)360 PARAMETER (cloia = 1.0E+02, cloib = -10.00, cloic = -0.6, cloid = 5.0) 366 361 ! cc PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.9, CLOID=5.0) 367 362 REAL rgammas 368 PARAMETER (rgammas =0.05)363 PARAMETER (rgammas = 0.05) 369 364 REAL crhl 370 PARAMETER (crhl =0.15)365 PARAMETER (crhl = 0.15) 371 366 ! cc PARAMETER (CRHL=0.70) 372 367 REAL t_coup 373 PARAMETER (t_coup =234.0)368 PARAMETER (t_coup = 234.0) 374 369 375 370 ! Variables locales: … … 392 387 END DO 393 388 394 DO k = 2, klev /2 - 1395 DO i = 1, klon 396 zdthdp = (t(i, k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) - &397 rd*0.5*(t(i,k)+t(i,k+1))/rcpd/paprs(i, k+1)398 zdthdp = zdthdp *cloia399 IF (pplay(i, k)>cetamb*paprs(i,1) .AND. zdthdp<zdthmin(i)) THEN389 DO k = 2, klev / 2 - 1 390 DO i = 1, klon 391 zdthdp = (t(i, k) - t(i, k + 1)) / (pplay(i, k) - pplay(i, k + 1)) - & 392 rd * 0.5 * (t(i, k) + t(i, k + 1)) / rcpd / paprs(i, k + 1) 393 zdthdp = zdthdp * cloia 394 IF (pplay(i, k)>cetamb * paprs(i, 1) .AND. zdthdp<zdthmin(i)) THEN 400 395 zdthmin(i) = zdthdp 401 396 invb(i) = k … … 407 402 kb = invb(i) 408 403 IF (thermcep) THEN 409 zdelta = max(0., sign(1., rtt-t(i,kb)))410 zqs = r2es *foeew(t(i,kb), zdelta)/pplay(i, kb)404 zdelta = max(0., sign(1., rtt - t(i, kb))) 405 zqs = r2es * foeew(t(i, kb), zdelta) / pplay(i, kb) 411 406 zqs = min(0.5, zqs) 412 zcor = 1. /(1.-retv*zqs)413 zqs = zqs *zcor407 zcor = 1. / (1. - retv * zqs) 408 zqs = zqs * zcor 414 409 ELSE 415 IF (t(i, kb)<t_coup) THEN416 zqs = qsats(t(i, kb))/pplay(i, kb)410 IF (t(i, kb)<t_coup) THEN 411 zqs = qsats(t(i, kb)) / pplay(i, kb) 417 412 ELSE 418 zqs = qsatl(t(i, kb))/pplay(i, kb)413 zqs = qsatl(t(i, kb)) / pplay(i, kb) 419 414 END IF 420 415 END IF 421 zcll = cloib*zdthmin(i) + cloic 422 zcll = min(1.0, max(0.0,zcll)) 423 zrhb = q(i, kb)/zqs 424 IF (zcll>0.0 .AND. zrhb<crhl) zcll = zcll*(1.-(crhl-zrhb)*cloid) 425 zcll = min(1.0, max(0.0,zcll)) 426 diafra(i, kb) = max(diafra(i,kb), zcll) 427 dialiq(i, kb) = diafra(i, kb)*rgammas*zqs 428 END DO 429 416 zcll = cloib * zdthmin(i) + cloic 417 zcll = min(1.0, max(0.0, zcll)) 418 zrhb = q(i, kb) / zqs 419 IF (zcll>0.0 .AND. zrhb<crhl) zcll = zcll * (1. - (crhl - zrhb) * cloid) 420 zcll = min(1.0, max(0.0, zcll)) 421 diafra(i, kb) = max(diafra(i, kb), zcll) 422 dialiq(i, kb) = diafra(i, kb) * rgammas * zqs 423 END DO 430 424 431 425 END SUBROUTINE diagcld2
Note: See TracChangeset
for help on using the changeset viewer.