Changeset 2984 for trunk/LMDZ.MARS/libf
- Timestamp:
- Jun 26, 2023, 5:44:59 PM (19 months ago)
- Location:
- trunk/LMDZ.MARS/libf/phymars
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/phymars/improvedclouds_mod.F
r2966 r2984 5 5 CONTAINS 6 6 7 subroutine improvedclouds(microtimestep, 8 & pplay,pt,pq, 9 & subpdqcloud,subpdtcloud, 10 & nq,tauscaling,mmean) 7 subroutine improvedclouds(ngrid,nlay,ptimestep, 8 & pplay,pt,pdt,pq,pdq,nq,tauscaling, 9 & imicro,zt,zq) 11 10 USE updaterad, ONLY: updaterice_micro, updaterccn 12 11 USE watersat_mod, ONLY: watersat … … 17 16 & igcm_hdo_vap,igcm_hdo_ice, 18 17 & qperemin 18 use conc_mod, only: mmean 19 19 use comcstfi_h, only: pi, cpp 20 20 use write_output_mod, only: write_output … … 41 41 c T. Navarro, debug,correction, new scheme (October-April 2011) 42 42 c A. Spiga, optimization (February 2012) 43 c J. Naar, from global to local (no more ngrid, nlay) to allow 44 c different microphysical timesteps (May 2023) 43 c J. Naar, adaptative subtimestep now done here (June 2023) 45 44 c------------------------------------------------------------------ 46 45 #include "callkeys.h" … … 49 48 c Inputs/outputs: 50 49 50 INTEGER, INTENT(IN) :: ngrid,nlay 51 51 INTEGER, INTENT(IN) :: nq ! nombre de traceurs 52 REAL, INTENT(IN) :: microtimestep! pas de temps physique (s)53 REAL, INTENT(IN) :: pplay ! pression au milieu des couches (Pa)54 REAL, INTENT(IN) :: pt ! temperature at the middle of the52 REAL, INTENT(IN) :: ptimestep ! pas de temps physique (s) 53 REAL, INTENT(IN) :: pplay(ngrid,nlay) ! pression au milieu des couches (Pa) 54 REAL, INTENT(IN) :: pt(ngrid,nlay) ! temperature at the middle of the 55 55 ! layers (K) 56 ! param. 57 REAL, INTENT(IN) :: pq(nq) ! traceur (kg/kg) 58 ! (kg/kg.s-1) 59 REAL, INTENT(IN) :: tauscaling ! Convertion factor for qdust and Ndust 60 REAL, INTENT(IN) :: mmean ! Mean atmospheric mass 61 62 REAL, INTENT(OUT) :: subpdqcloud(nq) ! tendance de la condensation 63 ! H2O(kg/kg.s-1) 64 REAL, INTENT(OUT) :: subpdtcloud ! tendance temperature due 65 ! a la chaleur latente 56 REAL, INTENT(IN) :: pdt(ngrid,nlay) ! temperature tendency (K/s) 57 REAL, INTENT(IN) :: pq(ngrid,nlay,nq) ! tracer (kg/kg) 58 REAL, INTENT(IN) :: pdq(ngrid,nlay,nq) ! tracer tendency (kg/kg/s) 59 REAL, INTENT(IN) :: tauscaling(ngrid) ! Convertion factor for qdust and Ndust 60 INTEGER, INTENT(IN) :: imicro ! nb. microphy calls(retrocompatibility) 61 62 REAL, INTENT(OUT) :: zq(ngrid,nlay,nq) ! tracers post microphy (kg/kg) 63 REAL, INTENT(OUT) :: zt(ngrid,nlay) ! temperature post microphy (K) 66 64 67 65 c------------------------------------------------------------------ … … 79 77 INTEGER ig,l,i 80 78 81 REAL zq(nq) ! local value of tracers 82 REAL zq0(nq) ! local initial value of tracers 83 REAL zt ! local value of temperature 84 REAL zqsat_tmp(1) ! saturation 85 REAL zqsat ! saturation 86 REAL lw !Latent heat of sublimation (J.kg-1) 79 REAL zqsat(ngrid,nlay) ! saturation 80 REAL lw !Latent heat of sublimation (J.kg-1) 87 81 REAL cste 88 82 REAL dMice ! mass of condensed ice 89 83 REAL dMice_hdo ! mass of condensed HDO ice 90 REAL alpha ! HDO equilibrium fractionation coefficient (Saturation=1)91 REAL alpha_c ! HDO real fractionation coefficient84 REAL alpha(ngrid,nlay) ! HDO equilibrium fractionation coefficient (Saturation=1) 85 REAL alpha_c(ngrid,nlay) ! HDO real fractionation coefficient 92 86 ! REAL sumcheck 93 87 REAL*8 ph2o ! Water vapor partial pressure (Pa) … … 105 99 REAL seq 106 100 107 REAL rice ! Ice mass mean radius (m)101 REAL rice(ngrid,nlay) ! Ice mass mean radius (m) 108 102 ! (r_c in montmessin_2004) 109 REAL rhocloud ! Cloud density (kg.m-3)110 REAL rdust ! Dust geometric mean radius (m)103 REAL rhocloud(ngrid,nlay) ! Cloud density (kg.m-3) 104 REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m) 111 105 112 106 REAL res ! Resistance growth … … 136 130 137 131 132 c---------------------------------- 133 c JN : used in subtimestep 134 REAL :: microtimestep! subdivision of ptimestep (s) 135 REAL :: subpdtcloud(ngrid,nlay) ! Temperature variation due to latent heat 136 REAL :: zq0(ngrid,nlay,nq) ! local initial value of tracers 137 c REAL zt0(ngrid,nlay,nq) ! local initial value of temperature 138 INTEGER :: zimicro(ngrid,nlay) ! Subdivision of ptimestep 139 INTEGER :: count_micro(ngrid,nlay) ! Number of microphys calls 140 REAL :: zpotcond_inst(ngrid,nlay) ! condensable water at the beginning of physics (kg/kg) 141 REAL :: zpotcond_full(ngrid,nlay) ! condensable water with integrated tendancies (kg/kg) 142 REAL :: zpotcond(ngrid,nlay) ! maximal condensable water (previous two) 143 REAL :: zqsat_tmp(1) ! maximal condensable water (previous two) 144 REAL :: spenttime ! time spent in while loop 145 REAL :: zdq ! used to compute adaptative timestep 146 LOGICAL :: ending_ts ! Condition to end while loop 147 138 148 139 149 c---------------------------------- … … 148 158 149 159 150 REAL satubf ,satuaf151 REAL res_out 160 REAL satubf(ngrid,nlay),satuaf(ngrid,nlay) 161 REAL res_out(ngrid,nlay) 152 162 153 163 … … 234 244 ! 1. Initialisation 235 245 !============================================================= 236 cste = 4*pi*rho_ice*microtimestep 237 238 res_out = 0 239 rice = 1.e-8 240 241 c Initialize the tendencies 242 subpdqcloud(1:nq)=0 243 subpdtcloud=0 244 245 c temperature and tracers previously incremented here 246 c now done outside, in watercloud_mod.F 247 c zt = pteff + sum_subpdt * microtimestep 248 c zq(1:nq) = pqeff(1:nq) + sum_subpdq(1:nq) * microtimestep 249 250 zq(1:nq)=pq(1:nq) 251 zt=pt 252 253 WHERE( zq(1:nq) < 1.e-30 ) 254 & zq(1:nq) = 1.e-30 255 256 zq0(1:nq) = zq(1:nq) 246 247 res_out(:,:) = 0 248 rice(:,:) = 1.e-8 249 250 c Initialize the temperature, tracers 251 zt(1:ngrid,1:nlay)=pt(1:ngrid,1:nlay) 252 zq(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq) 253 subpdtcloud(1:ngrid,1:nlay)=0 254 255 WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 ) 256 & zq(1:ngrid,1:nlay,1:nq) = 1.e-30 257 257 258 258 259 !============================================================= … … 263 264 dev2 = 1. / ( sqrt(2.) * sigma_ice ) 264 265 265 call watersat(1,(/zt/),(/pplay/),zqsat_tmp) 266 zqsat=zqsat_tmp(1) 266 c Compute the condensable amount of water vapor 267 c or the sublimable water ice (if negative) 268 call watersat(ngrid*nlay,zt+pdt*ptimestep,pplay,zqsat) 269 zpotcond_full=(zq(:,:,igcm_h2o_vap)+ 270 & pdq(:,:,igcm_h2o_vap)*ptimestep) - zqsat 271 zpotcond_full=max(zpotcond_full,-zq(:,:,igcm_h2o_ice)) 272 call watersat(ngrid*nlay,zt,pplay,zqsat) 273 zpotcond_inst=zq(:,:,igcm_h2o_vap) - zqsat 274 zpotcond_inst=max(zpotcond_inst,-zq(:,:,igcm_h2o_ice)) 275 c zpotcond is the most strict criterion between the two 276 DO l=1,nlay 277 DO ig=1,ngrid 278 if (zpotcond_full(ig,l).gt.0.) then 279 zpotcond(ig,l)=max(zpotcond_inst(ig,l),zpotcond_full(ig,l)) 280 else if (zpotcond_full(ig,l).le.0.) then 281 zpotcond(ig,l)=min(zpotcond_inst(ig,l),zpotcond_full(ig,l)) 282 endif! (zpotcond_full.gt.0.) then 283 ENDDO !ig=1,ngrid 284 ENDDO !l=1,nlay 285 267 286 countcells = 0 268 287 zimicro(1:ngrid,1:nlay)=imicro 288 count_micro(1:ngrid,1:nlay)=0 289 290 c Main loop over the GCM's grid 291 DO l=1,nlay 292 DO ig=1,ngrid 293 c Subtimestep : here we go 294 IF (cloud_adapt_ts) then 295 call adapt_imicro(ptimestep,zpotcond(ig,l), 296 & zimicro(ig,l)) 297 ENDIF! (cloud_adapt_ts) then 298 spenttime = 0. 299 ending_ts=.false. 300 print*, "ig,l :", ig, l 301 DO while (.not.ending_ts) 302 microtimestep=ptimestep/real(zimicro(ig,l)) 303 zq0(1:ngrid,1:nlay,1:nq) = zq 304 305 ! Check if we are integrating over ptimestep 306 if (spenttime+microtimestep.ge.ptimestep) then 307 microtimestep=ptimestep-spenttime 308 ! If so : last call ! 309 ending_ts=.true. 310 endif! (spenttime+microtimestep.ge.ptimestep) then 311 call watersat(1,(/zt(ig,l)/),(/pplay(ig,l)/),zqsat_tmp) 312 zqsat(ig,l)=zqsat_tmp(1) 269 313 c Get the partial pressure of water vapor and its saturation ratio 270 ph2o = zq(igcm_h2o_vap) * (mmean/18.) * pplay271 satu = zq(igcm_h2o_vap) / zqsat314 ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l) 315 satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l) 272 316 273 317 !============================================================= … … 275 319 !============================================================= 276 320 277 IF ( satu .ge. 1. ) THEN ! if there is condensation278 279 call updaterccn(zq(igcm_dust_mass),280 & zq(igcm_dust_number),rdust,tauscaling)321 IF ( satu .ge. 1. ) THEN ! if there is condensation 322 323 call updaterccn(zq(ig,l,igcm_dust_mass), 324 & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig)) 281 325 282 326 283 327 c Expand the dust moments into a binned distribution 284 Mo = zq(ig cm_dust_mass)* tauscaling+ 1.e-30285 No = zq(ig cm_dust_number)* tauscaling+ 1.e-30286 Rn = rdust 328 Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) + 1.e-30 329 No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30 330 Rn = rdust(ig,l) 287 331 Rn = -log(Rn) 288 332 Rm = Rn - 3. * sigma_ice*sigma_ice … … 324 368 325 369 c Get the rates of nucleation 326 call nuclea(ph2o,zt ,satu,n_aer,rate)370 call nuclea(ph2o,zt(ig,l),satu,n_aer,rate) 327 371 328 372 dN = 0. … … 335 379 336 380 c Update Dust particles 337 zq(ig cm_dust_mass) =338 & zq(ig cm_dust_mass) + dM/ tauscaling !max(tauscaling,1.e-10)339 zq(ig cm_dust_number) =340 & zq(ig cm_dust_number) + dN/ tauscaling !max(tauscaling,1.e-10)381 zq(ig,l,igcm_dust_mass) = 382 & zq(ig,l,igcm_dust_mass) + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10) 383 zq(ig,l,igcm_dust_number) = 384 & zq(ig,l,igcm_dust_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10) 341 385 c Update CCNs 342 zq(ig cm_ccn_mass) =343 & zq(ig cm_ccn_mass) - dM/ tauscaling !max(tauscaling,1.e-10)344 zq(ig cm_ccn_number) =345 & zq(ig cm_ccn_number) - dN/ tauscaling !max(tauscaling,1.e-10)386 zq(ig,l,igcm_ccn_mass) = 387 & zq(ig,l,igcm_ccn_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10) 388 zq(ig,l,igcm_ccn_number) = 389 & zq(ig,l,igcm_ccn_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10) 346 390 347 ENDIF ! of is satu >1391 ENDIF ! of is satu >1 348 392 349 393 !============================================================= … … 355 399 c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible. 356 400 357 IF ( zq(ig cm_ccn_number)*tauscaling.ge. 1.) THEN ! we trigger crystal growth401 IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth 358 402 359 403 360 call updaterice_micro(zq(ig cm_h2o_ice),361 & zq(ig cm_ccn_mass),zq(igcm_ccn_number),362 & tauscaling ,rice,rhocloud)363 364 No = zq(ig cm_ccn_number)* tauscaling+ 1.e-30404 call updaterice_micro(zq(ig,l,igcm_h2o_ice), 405 & zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number), 406 & tauscaling(ig),rice(ig,l),rhocloud(ig,l)) 407 408 No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30 365 409 366 410 c saturation at equilibrium 367 411 c rice should not be too small, otherwise seq value is not valid 368 seq = exp(2.*sig(zt )*mh2o / (rho_ice*rgp*zt*369 & max(rice ,1.e-7)))412 seq = exp(2.*sig(zt(ig,l))*mh2o / (rho_ice*rgp*zt(ig,l)* 413 & max(rice(ig,l),1.e-7))) 370 414 371 415 c get resistance growth 372 call growthrate(zt ,pplay,373 & real(ph2o/satu),rice ,res,Dv)374 375 res_out = res416 call growthrate(zt(ig,l),pplay(ig,l), 417 & real(ph2o/satu),rice(ig,l),res,Dv) 418 419 res_out(ig,l) = res 376 420 377 421 ccccccc implicit scheme of mass growth 422 c cste here must be computed at each step 423 cste = 4*pi*rho_ice*microtimestep 424 c print*, 'ig',ig 425 c print*, 'l',l 426 c print*, 'cste',cste 427 c print*, 'seq',seq 428 c print*, 'zqsat',zqsat(ig,l) 429 c print*, 'zq vap',zq(ig,l,igcm_h2o_vap) 430 c print*, 'res',res 431 c print*, 'No',No 432 c print*, 'rice',rice(ig,l) 378 433 379 434 dMice = 380 & (zq(ig cm_h2o_vap)-seq*zqsat)381 & /(res*zqsat /(cste*No*rice) + 1.)435 & (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l)) 436 & /(res*zqsat(ig,l)/(cste*No*rice(ig,l)) + 1.) 382 437 383 438 384 439 ! With the above scheme, dMice cannot be bigger than vapor, 385 440 ! but can be bigger than all available ice. 386 dMice = max(dMice,-zq(ig cm_h2o_ice))387 dMice = min(dMice,zq(ig cm_h2o_vap)) ! this should be useless...388 389 zq(ig cm_h2o_ice) = zq(igcm_h2o_ice)+dMice390 zq(ig cm_h2o_vap) = zq(igcm_h2o_vap)-dMice441 dMice = max(dMice,-zq(ig,l,igcm_h2o_ice)) 442 dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless... 443 444 zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice 445 zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice 391 446 392 447 … … 394 449 395 450 ! latent heat release 396 lw=(2834.3-0.28*(zt -To)-397 & 0.004*(zt -To)*(zt-To))*1.e+3398 subpdtcloud = dMice*lw/cpp/microtimestep451 lw=(2834.3-0.28*(zt(ig,l)-To)- 452 & 0.004*(zt(ig,l)-To)*(zt(ig,l)-To))*1.e+3 453 subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep 399 454 400 455 c Special case of the isotope of water HDO … … 405 460 if (hdofrac) then 406 461 !! Calculation of the HDO vapor coefficient 407 Dv_hdo = 1./3. * sqrt( 8*kbz*zt /(pi*mhdo/nav) )408 & * kbz * zt /409 & ( pi * pplay * (molco2+molhdo)*(molco2+molhdo)462 Dv_hdo = 1./3. * sqrt( 8*kbz*zt(ig,l)/(pi*mhdo/nav) ) 463 & * kbz * zt(ig,l) / 464 & ( pi * pplay(ig,l) * (molco2+molhdo)*(molco2+molhdo) 410 465 & * sqrt(1.+mhdo/mco2) ) 411 466 !! Calculation of the fractionnation coefficient at equilibrium 412 c alpha = exp(16288./zt**2.-9.34d-2) ! Merlivat and Nief et al. 1967 413 alpha = exp(13525./zt**2.-5.59d-2) ! Lamb et al. 2017 467 alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2) 468 c alpha = exp(13525./zt(ig,l)**2.-5.59d-2) !Lamb 414 469 !! Calculation of the 'real' fractionnation coefficient 415 alpha_c = (alpha*satu)/416 & ( (alpha *(Dv/Dv_hdo)*(satu-1.)) + 1.)417 c alpha_c = alpha! to test without the effect of cinetics470 alpha_c(ig,l) = (alpha(ig,l)*satu)/ 471 & ( (alpha(ig,l)*(Dv/Dv_hdo)*(satu-1.)) + 1.) 472 c alpha_c(ig,l) = alpha(ig,l) ! to test without the effect of cinetics 418 473 else 419 alpha_c = 1.d0474 alpha_c(ig,l) = 1.d0 420 475 endif 421 if (zq0(ig cm_h2o_vap).gt.qperemin) then476 if (zq0(ig,l,igcm_h2o_vap).gt.qperemin) then 422 477 dMice_hdo= 423 & dMice*alpha_c *424 & ( zq0(ig cm_hdo_vap)425 & /zq0(ig cm_h2o_vap) )478 & dMice*alpha_c(ig,l)* 479 & ( zq0(ig,l,igcm_hdo_vap) 480 & /zq0(ig,l,igcm_h2o_vap) ) 426 481 else 427 482 dMice_hdo=0. … … 429 484 !! sublimation 430 485 else 431 if (zq0(ig cm_h2o_ice).gt.qperemin) then486 if (zq0(ig,l,igcm_h2o_ice).gt.qperemin) then 432 487 dMice_hdo= 433 488 & dMice* 434 & ( zq0(ig cm_hdo_ice)435 & /zq0(ig cm_h2o_ice) )489 & ( zq0(ig,l,igcm_hdo_ice) 490 & /zq0(ig,l,igcm_h2o_ice) ) 436 491 else 437 492 dMice_hdo=0. … … 439 494 endif !if (dMice.gt.0.0) 440 495 441 dMice_hdo = max(dMice_hdo,-zq(ig cm_hdo_ice))442 dMice_hdo = min(dMice_hdo,zq(ig cm_hdo_vap))443 444 zq(ig cm_hdo_ice) = zq(igcm_hdo_ice)+dMice_hdo445 zq(ig cm_hdo_vap) = zq(igcm_hdo_vap)-dMice_hdo496 dMice_hdo = max(dMice_hdo,-zq(ig,l,igcm_hdo_ice)) 497 dMice_hdo = min(dMice_hdo,zq(ig,l,igcm_hdo_vap)) 498 499 zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+dMice_hdo 500 zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)-dMice_hdo 446 501 447 502 endif ! if (hdo) … … 453 508 c If all the ice particles sublimate, all the condensation 454 509 c nuclei are released: 455 if (zq(ig cm_h2o_ice).le.1.e-28) then510 if (zq(ig,l,igcm_h2o_ice).le.1.e-28) then 456 511 457 512 c Water 458 zq(ig cm_h2o_vap) = zq(igcm_h2o_vap)459 & + zq(ig cm_h2o_ice)460 zq(ig cm_h2o_ice) = 0.513 zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap) 514 & + zq(ig,l,igcm_h2o_ice) 515 zq(ig,l,igcm_h2o_ice) = 0. 461 516 if (hdo) then 462 zq(ig cm_hdo_vap) = zq(igcm_hdo_vap)463 & + zq(ig cm_hdo_ice)464 zq(ig cm_hdo_ice) = 0.517 zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap) 518 & + zq(ig,l,igcm_hdo_ice) 519 zq(ig,l,igcm_hdo_ice) = 0. 465 520 endif 466 521 c Dust particles 467 zq(ig cm_dust_mass) = zq(igcm_dust_mass)468 & + zq(ig cm_ccn_mass)469 zq(ig cm_dust_number) = zq(igcm_dust_number)470 & + zq(ig cm_ccn_number)522 zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) 523 & + zq(ig,l,igcm_ccn_mass) 524 zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) 525 & + zq(ig,l,igcm_ccn_number) 471 526 c CCNs 472 zq(ig cm_ccn_mass) = 0.473 zq(ig cm_ccn_number) = 0.527 zq(ig,l,igcm_ccn_mass) = 0. 528 zq(ig,l,igcm_ccn_number) = 0. 474 529 475 530 endif 476 531 477 532 ENDIF !of if Nccn>1 533 478 534 479 480 ! Get cloud tendencies 481 subpdqcloud(igcm_h2o_vap) = 482 & (zq(igcm_h2o_vap) - 483 & zq0(igcm_h2o_vap))/microtimestep 484 subpdqcloud(igcm_h2o_ice) = 485 & (zq(igcm_h2o_ice) - 486 & zq0(igcm_h2o_ice))/microtimestep 487 if (hdo) then 488 subpdqcloud(igcm_hdo_vap) = 489 & (zq(igcm_hdo_vap) - 490 & zq0(igcm_hdo_vap))/microtimestep 491 subpdqcloud(igcm_hdo_ice) = 492 & (zq(igcm_hdo_ice) - 493 & zq0(igcm_hdo_ice))/microtimestep 494 endif 495 subpdqcloud(igcm_ccn_mass) = 496 & (zq(igcm_ccn_mass) - 497 & zq0(igcm_ccn_mass))/microtimestep 498 subpdqcloud(igcm_ccn_number) = 499 & (zq(igcm_ccn_number) - 500 & zq0(igcm_ccn_number))/microtimestep 501 502 if (scavenging) then 503 504 subpdqcloud(igcm_dust_mass) = 505 & (zq(igcm_dust_mass) - 506 & zq0(igcm_dust_mass))/microtimestep 507 subpdqcloud(igcm_dust_number) = 508 & (zq(igcm_dust_number) - 509 & zq0(igcm_dust_number))/microtimestep 510 511 endif 535 ! No more getting tendency : we increment tracers & temp directly 536 537 ! Increment tracers & temp for following microtimestep 538 ! Dust : 539 ! Special treatment of dust_mass / number for scavenging ? 540 IF (scavenging) THEN 541 zq(ig,l,igcm_dust_mass) =zq(ig,l,igcm_dust_mass)+ 542 & pdq(ig,l,igcm_dust_mass)*microtimestep 543 zq(ig,l,igcm_dust_number) =zq(ig,l,igcm_dust_number)+ 544 & pdq(ig,l,igcm_dust_number)*microtimestep 545 ELSE 546 zq(ig,l,igcm_dust_mass) =zq0(ig,l,igcm_dust_mass) 547 zq(ig,l,igcm_dust_number) =zq0(ig,l,igcm_dust_number) 548 ENDIF !(scavenging) THEN 549 zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + 550 & pdq(ig,l,igcm_ccn_mass)*microtimestep 551 zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + 552 & pdq(ig,l,igcm_ccn_number)*microtimestep 553 554 ! Water : 555 zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+ 556 & pdq(ig,l,igcm_h2o_ice)*microtimestep 557 zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)+ 558 & pdq(ig,l,igcm_h2o_vap)*microtimestep 559 560 ! HDO (if computed) : 561 IF (hdo) THEN 562 zq(ig,l,igcm_hdo_ice) = 563 & zq(ig,l,igcm_hdo_ice)+ 564 & pdq(ig,l,igcm_hdo_ice)*microtimestep 565 zq(ig,l,igcm_hdo_vap) = 566 & zq(ig,l,igcm_hdo_vap)+ 567 & pdq(ig,l,igcm_hdo_vap)*microtimestep 568 ENDIF ! hdo 569 570 c Could also set subpdtcloud to 0 if not activice to make it simpler 571 c or change name of the flag 572 IF (.not.activice) THEN 573 subpdtcloud(ig,l)=0. 574 ENDIF 575 ! Temperature : 576 zt(ig,l) = zt(ig,l)+(pdt(ig,l)+ 577 & subpdtcloud(ig,l))*microtimestep 578 579 c Prevent negative tracers 580 WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 ) 581 & zq(1:ngrid,1:nlay,1:nq) = 1.e-30 582 583 c Increment time spent in here 584 spenttime=spenttime+microtimestep 585 count_micro(ig,l)=count_micro(ig,l)+1 586 IF (cloud_adapt_ts) then 587 c Estimation of how much is actually condensing/subliming 588 c print*, 'zq',zq(ig,l,igcm_h2o_vap) 589 c print*, 'zq0',zq0(ig,l,igcm_h2o_vap) 590 c print*, 'ptimestep',ptimestep 591 zdq=(zq(ig,l,igcm_h2o_vap)-zq0(ig,l,igcm_h2o_vap)) 592 & *ptimestep/microtimestep 593 c Estimation of how much is actually condensing/subliming 594 c zdq=min(abs(dMice)*ptimestep,zpotcond(ig,l)) 595 call adapt_imicro(ptimestep,zdq, 596 & zimicro(ig,l)) 597 ENDIF! (cloud_adapt_ts) then 598 ENDDO ! while (.not. ending_ts) 599 ENDDO ! of ig loop 600 ENDDO ! of nlayer loop 601 602 603 c------ Useful outputs to check how it went 604 call write_output("zpotcond_inst","zpotcond_inst "// 605 & "microphysics","(kg/kg)",zpotcond_inst(:,:)) 606 call write_output("zpotcond_full","zpotcond_full "// 607 & "microphysics","(kg/kg)",zpotcond_full(:,:)) 608 call write_output("zpotcond","zpotcond "// 609 & "microphysics","(kg/kg)",zpotcond(:,:)) 610 call write_output("count_micro","count_micro "// 611 & "after microphysics","integer",count_micro(:,:)) 512 612 513 613 … … 519 619 520 620 ! error2d(:) = 0. 521 ! error2d(ig) = max(abs(error_out),error2d(ig)) 522 satubf = zq0(igcm_h2o_vap)/zqsat 523 satuaf = zq(igcm_h2o_vap)/zqsat 524 525 c print*, 'count is ',countcells, ' i.e. ', 526 c & countcells*100/(nlay*ngrid), '% for microphys computation' 621 DO l=1,nlay 622 DO ig=1,ngrid 623 ! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig)) 624 satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l) 625 satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l) 626 ENDDO 627 ENDDO 628 629 print*, 'count is ',countcells, ' i.e. ', 630 & countcells*100/(nlay*ngrid), '% for microphys computation' 527 631 528 632 #ifndef MESOSCALE 529 633 ! IF (ngrid.ne.1) THEN ! 3D 530 ! call write_output("satu","ratio saturation","",634 ! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3, 531 635 ! & satu_out) 532 ! call write_output("dM","ccn variation","kg/kg",636 ! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3, 533 637 ! & dM_out) 534 ! call write_output("dN","ccn variation","#",638 ! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3, 535 639 ! & dN_out) 536 ! call write_output("error","dichotomy max error","%",640 ! call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2, 537 641 ! & error2d) 538 ! call write_output("zqsat","zqsat","kg",642 ! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3, 539 643 ! & zqsat) 540 644 ! ENDIF 541 645 542 646 ! IF (ngrid.eq.1) THEN ! 1D 543 ! call write_output("error","incertitude sur glace","%",647 ! call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1, 544 648 ! & error_out) 545 call write_output("resist","resistance","s/m2",649 call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1, 546 650 & res_out) 547 call write_output("satu_bf","satu before","kg/kg",651 call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1, 548 652 & satubf) 549 call write_output("satu_af","satu after","kg/kg",653 call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1, 550 654 & satuaf) 551 call write_output("vapbf","h2ovap before","kg/kg",552 & zq0( igcm_h2o_vap))553 call write_output("vapaf","h2ovap after","kg/kg",554 & zq( igcm_h2o_vap))555 call write_output("icebf","h2oice before","kg/kg",556 & zq0( igcm_h2o_ice))557 call write_output("iceaf","h2oice after","kg/kg",558 & zq( igcm_h2o_ice))559 call write_output("ccnbf","ccn before","/kg",560 & zq0( igcm_ccn_number))561 call write_output("ccnaf","ccn after","/kg",562 & zq( igcm_ccn_number))563 c call write_output("growthrate","growth rate","m^2/s",655 call WRITEdiagfi(ngrid,"vapbf","h2ovap before","kg/kg",1, 656 & zq0(1,1,igcm_h2o_vap)) 657 call WRITEdiagfi(ngrid,"vapaf","h2ovap after","kg/kg",1, 658 & zq(1,1,igcm_h2o_vap)) 659 call WRITEdiagfi(ngrid,"icebf","h2oice before","kg/kg",1, 660 & zq0(1,1,igcm_h2o_ice)) 661 call WRITEdiagfi(ngrid,"iceaf","h2oice after","kg/kg",1, 662 & zq(1,1,igcm_h2o_ice)) 663 call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1, 664 & zq0(1,1,igcm_ccn_number)) 665 call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1, 666 & zq(1,1,igcm_ccn_number)) 667 c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1, 564 668 c & gr_out) 565 c call write_output("nuclearate","nucleation rate","",669 c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1, 566 670 c & rate_out) 567 c call write_output("dM","ccn variation","kg",671 c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1, 568 672 c & dM_out) 569 c call write_output("dN","ccn variation","#",673 c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1, 570 674 c & dN_out) 571 call write_output("zqsat","p vap sat","kg/kg",675 call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1, 572 676 & zqsat) 573 ! call write_output("satu","ratio saturation","",574 ! & satu_out (:,:))575 call write_output("rice","ice radius","m",677 ! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1, 678 ! & satu_out) 679 call WRITEdiagfi(ngrid,"rice","ice radius","m",1, 576 680 & rice) 577 ! call write_output("rdust_sca","rdust","m",681 ! call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1, 578 682 ! & rdust) 579 ! call write_output("rsedcloud","rsedcloud","m",683 ! call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1, 580 684 ! & rsedcloud) 581 ! call write_output("rhocloud","rhocloud","kg.m-3",685 ! call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1, 582 686 ! & rhocloud) 583 687 ! ENDIF … … 636 740 637 741 END SUBROUTINE improvedclouds 742 743 SUBROUTINE adapt_imicro(ptimestep,potcond, 744 $ zimicro) 745 746 c Adaptative timestep for water ice clouds. 747 c Works using a powerlaw to compute the minimal duration of subtimestep 748 c (in s) should all the avalaible vapor (resp. ice) condenses (resp.sublimates) 749 c Then, we use the instantaneous vap (ice) gradient extrapolated to the 750 c rest of duration to increase subtimestep duration, for computing 751 c efficiency. 752 753 real,intent(in) :: ptimestep ! total duration of physics (sec) 754 real,intent(in) :: potcond ! condensible vapor / sublimable ice(kg/kg) 755 real :: alpha, beta ! Coefficients for power law 756 integer,intent(out) :: zimicro ! number of ptimestep division 757 758 c zimicro = 30 759 c Coefficients good enough for present-day Mars : 760 alpha = 1.87485684e+09 761 beta = 1.45655856e+00 762 c Coefficients covering high obliquity scenarios : 763 c alpha=3.36711332e+15 764 c beta=1.98636414e+00 765 c nimicro=min(max(alpha*abs(potcond)**beta,5.),7000.) 766 c zimicro=ceiling((ptimestep/timeleft)*nimicro) 767 c zimicro=ceiling((timeleft/ptimestep)*nimicro) 768 zimicro=ceiling(min(max(alpha*abs(potcond)**beta,5.),7000.)) 769 770 END SUBROUTINE adapt_imicro 638 771 639 772 END MODULE improvedclouds_mod -
trunk/LMDZ.MARS/libf/phymars/watercloud_mod.F
r2966 r2984 47 47 c 2004 - 2012 48 48 c 49 c 2023: J. Naar, adding different subtimestep on each grid cell 50 c plus not doing microphysics if no water present 51 c plus simpleclouds no longer in the loop for microphysics 49 c 2023: J. Naar, now with adaptative timestep for improvedclouds 50 c (done in improvedclouds_mod). 52 51 c======================================================================= 53 52 … … 146 145 ! Scheme for adaptative timestep J. Naar 2023 147 146 c LOGICAL :: computed_micro(ngrid,nlay) ! Check if microphy was done in this cell 148 REAL :: co mputed_micro(ngrid,nlay) ! Check if microphy was done inthis cell (logical)147 REAL :: count_micro(ngrid,nlay) ! Initially computed microtimestep 149 148 REAL :: zt_micro(ngrid,nlay) ! Temperature during microphysics (K) 150 149 REAL :: zq_micro(ngrid,nlay,nq) ! Tracers during microphysics (kg/kg) … … 155 154 REAL :: zpotcond(ngrid,nlay) ! maximal condensable water, used to 156 155 compute adaptative subdivision of ptimestep 156 REAL :: spenttime ! timespent 157 REAL :: zdq ! used to compute adaptative timestep 157 158 158 159 … … 298 299 END IF ! end if (CLFvarying) 299 300 c------------------------------------------------------------------ 300 c Time subsampling for microphysics301 c Cloud physics (nucleation, condensation / sublimation) 301 302 c------------------------------------------------------------------ 302 303 rhocloud(1:ngrid,1:nlay) = rho_dust 303 304 304 305 c Initialisation of all the stuff JN2023 306 c computed_micro(1:ngrid,1:nlay)=.false. 307 computed_micro(1:ngrid,1:nlay)=0. 305 c Initialisation of all the stuff (JN,2023) 308 306 zt_micro(:,:)=pt(:,:) 309 307 zq_micro(:,:,:)=pq(:,:,:) 310 zq_micro(:,:,:)=pq(:,:,:) 311 call watersat(ngrid*nlay,zt_micro,pplay,zqsat_micro) 312 zpotcond_inst=zq_micro(:,:,igcm_h2o_vap) - zqsat_micro 313 call watersat(ngrid*nlay,zt_micro+pdt*ptimestep,pplay,zqsat_micro) 314 zpotcond_full=(zq_micro(:,:,igcm_h2o_vap)+ 315 & pdq(:,:,igcm_h2o_vap)*ptimestep) - zqsat_micro 316 zimicro(1:ngrid,1:nlay)=imicro 317 if (cloud_adapt_ts) then 318 call adapt_imicro(ptimestep,zpotcond(ig,l), 319 $ zimicro(ig,l)) 320 endif! (cloud_adapt_ts) then 321 DO l=1,nlay 322 DO ig=1,ngrid 323 c Start by computing the condensable water vapor amount 324 if (zpotcond_full(ig,l).gt.0.) then 325 zpotcond(ig,l)=max(zpotcond_inst(ig,l),zpotcond_full(ig,l)) 326 else if (zpotcond_full(ig,l).le.0.) then 327 zpotcond(ig,l)=min(zpotcond_inst(ig,l),zpotcond_full(ig,l)) 328 endif! (zpotcond_full.gt.0.) then 329 microtimestep=ptimestep/real(zimicro(ig,l)) 330 c Check if microphysics is even needed, that is if enough action 331 c is happening water-wise 332 if ((pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep 333 & .gt.1e-22) .or. (abs(zpotcond(ig,l)).gt.1e-22)) then 334 c Eventuellement sortir simpleclouds de la boucle egalement 335 computed_micro(ig,l)=1. 336 DO microstep=1,zimicro(ig,l) 337 338 339 c JN : incrementing after main microphysics scheme 340 c Previously we were incrementing tendencies, we now 341 c increment tracers and temperature directly 342 c We are thus starting at the end of the first iteration 343 c 308 344 309 c------------------------------------------------------------------- 345 310 c 1. Main call to the different cloud schemes: 346 311 c------------------------------------------------ 347 IF (microphys) THEN 348 CALL improvedclouds(microtimestep, 349 & pplay(ig,l),zt_micro(ig,l), 350 & zq_micro(ig,l,:),subpdqcloud(ig,l,:), 351 & subpdtcloud(ig,l),nq,tauscaling(ig),mmean(ig,l)) 352 353 ELSE 354 c Simpleclouds should maybe be taken out and put in a specific loop ? 355 CALL simpleclouds(ngrid,nlay,microtimestep, 312 c ds. 313 IF (microphys) THEN 314 CALL improvedclouds(ngrid,nlay,ptimestep, 315 & pplay,pt,pdt,pq,pdq,nq,tauscaling,imicro, 316 & zt_micro,zq_micro) 317 318 ELSE 319 320 c Specific loop for simpleclouds. 321 DO l=1,nlay 322 DO ig=1,ngrid 323 CALL simpleclouds(ngrid,nlay,ptimestep, 356 324 & pplay,pzlay,pteff,sum_subpdt, 357 325 & pqeff,sum_subpdq,subpdqcloud,subpdtcloud, 358 326 & nq,tau,rice) 359 ENDIF360 361 327 c------------------------------------------------------------------- 362 328 c 2. Updating tracers and temperature after cloud scheme: 329 c For improved clouds (with microphysics) this is done directly 330 c in the microphysics, during the subtimestep 331 c I put it like that to be retrocompatible (JN) 363 332 c----------------------------------------------- 364 333 365 IF (microphys) THEN366 zq_micro(ig,l,igcm_dust_mass) =367 & zq_micro(ig,l,igcm_dust_mass)+(pdq(ig,l,igcm_dust_mass)368 & +subpdqcloud(ig,l,igcm_dust_mass))*microtimestep369 zq_micro(ig,l,igcm_dust_number) =370 & zq_micro(ig,l,igcm_dust_number)371 & +(pdq(ig,l,igcm_dust_number)372 & + subpdqcloud(ig,l,igcm_dust_number))*microtimestep373 zq_micro(ig,l,igcm_ccn_mass) =374 & zq_micro(ig,l,igcm_ccn_mass) +375 & (pdq(ig,l,igcm_ccn_mass)376 & +subpdqcloud(ig,l,igcm_ccn_mass))*microtimestep377 zq_micro(ig,l,igcm_ccn_number) =378 & zq_micro(ig,l,igcm_ccn_number) +379 & (pdq(ig,l,igcm_ccn_number)380 & + subpdqcloud(ig,l,igcm_ccn_number))*microtimestep381 ENDIF382 334 zq_micro(ig,l,igcm_h2o_ice) = 383 335 & zq_micro(ig,l,igcm_h2o_ice)+ 384 336 & (pdq(ig,l,igcm_h2o_ice) 385 & + subpdqcloud(ig,l,igcm_h2o_ice))* microtimestep337 & + subpdqcloud(ig,l,igcm_h2o_ice))*ptimestep 386 338 zq_micro(ig,l,igcm_h2o_vap) = 387 339 & zq_micro(ig,l,igcm_h2o_vap)+ 388 340 & (pdq(ig,l,igcm_h2o_vap) 389 & + subpdqcloud(ig,l,igcm_h2o_vap))* microtimestep341 & + subpdqcloud(ig,l,igcm_h2o_vap))*ptimestep 390 342 391 343 IF (hdo) THEN … … 393 345 & zq_micro(ig,l,igcm_hdo_ice)+ 394 346 & (pdq(ig,l,igcm_hdo_ice) 395 & + subpdqcloud(ig,l,igcm_hdo_ice))* microtimestep347 & + subpdqcloud(ig,l,igcm_hdo_ice))*ptimestep 396 348 zq_micro(ig,l,igcm_hdo_vap) = 397 349 & zq_micro(ig,l,igcm_hdo_vap)+ 398 350 & (pdq(ig,l,igcm_hdo_vap) 399 & + subpdqcloud(ig,l,igcm_hdo_vap))* microtimestep351 & + subpdqcloud(ig,l,igcm_hdo_vap))*ptimestep 400 352 ENDIF ! hdo 401 353 402 354 c Could also set subpdtcloud to 0 if not activice to make it simpler 403 zt_micro(ig,l) = zt_micro(ig,l)+ 404 & pdt(ig,l)*microtimestep 405 IF (activice) THEN 406 zt_micro(ig,l) = zt_micro(ig,l)+ 407 & subpdtcloud(ig,l)*microtimestep 408 ENDIF 409 ! !! Example of how to use writediagmicrofi useful to 410 ! !! get outputs at each microphysical sub-timestep (better to be used in 1D) 411 ! CALL WRITEDIAGMICROFI(ngrid,imicro,microstep, 412 ! & microtimestep,'subpdtcloud', 413 ! & 'subpdtcloud','K/s',1,subpdtcloud(:,:)) 414 415 ENDDO ! of DO microstep=1,imicro 416 endif! (zq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep 417 ! & .gt.1e-22).or.(abs(zpotcond).gt.1e-22) then 418 ENDDO ! ig=1,ngrid 419 ENDDO ! l=1,nlay 420 421 422 c------ Useful outputs to check how it went 423 call write_output("computed_micro","computed_micro "// 424 & "after microphysics","logical",computed_micro(:,:)) 425 call write_output("zimicro","Used number of subtimestep "// 426 & "in cloud microphysics","integer",real(zimicro(:,:))) 355 c or change name of the flag 356 IF (activice) THEN 357 zt_micro(ig,l) = zt_micro(ig,l)+ 358 & subpdtcloud(ig,l)*ptimestep 359 ENDIF 360 361 ENDDO !ig=1,ngrid 362 ENDDO !l=1,nlay 363 ENDIF 364 365 366 427 367 c------------------------------------------------------------------- 428 368 c 3. Compute final tendencies after time loop: … … 786 726 787 727 788 SUBROUTINE adapt_imicro(ptimestep,potcond,789 $ zimicro)790 791 c Pas de temps adaptatif pour les nuages792 793 real,intent(in) :: ptimestep ! total duration of physics (sec)794 real,intent(in) :: potcond ! total duration of physics (sec)795 real :: alpha, beta ! total duration of physics (sec)796 integer,intent(out) :: zimicro ! number of ptimestep division797 798 c zimicro = ptimestep*alpha*potcond**beta799 zimicro = 30800 801 END SUBROUTINE adapt_imicro802 803 804 728 END MODULE watercloud_mod
Note: See TracChangeset
for help on using the changeset viewer.