Changeset 522
- Timestamp:
- Feb 10, 2012, 12:48:13 PM (13 years ago)
- File:
-
- 1 edited
-
trunk/LMDZ.MARS/libf/phymars/watercloud.F (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/phymars/watercloud.F
r459 r522 1 SUBROUTINE watercloud(ngrid,nlay, 1 SUBROUTINE watercloud(ngrid,nlay,ptimestep, 2 2 & pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt, 3 3 & pq,pdq,pdqcloud,pdqscloud,pdtcloud, 4 4 & nq,tau,tauscaling,rdust,rice,nuice, 5 5 & rsedcloud,rhocloud) 6 ! to use 'getin' 7 USE ioipsl_getincom 6 8 IMPLICIT NONE 7 9 … … 13 15 c - An improved microphysical scheme (see improvedclouds.F) 14 16 c 17 c There is a time loop specific to cloud formation and sedimentation 18 c due to timescales smaller than the GCM integration timestep. 19 c 15 20 c Authors: Franck Montmessin, Francois Forget, Ehouarn Millour, 16 c J.-B. Madeleine 21 c J.-B. Madeleine, Thomas Navarro 17 22 c 18 c 2004 - Oct. 201123 c 2004 - 2012 19 24 c======================================================================= 20 25 … … 35 40 36 41 INTEGER ngrid,nlay 37 integernq ! nombre de traceurs42 INTEGER nq ! nombre de traceurs 38 43 REAL ptimestep ! pas de temps physique (s) 39 44 REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa) 40 45 REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa) 41 REAL pdpsrf(ngrid) ! tend ance surf pressure46 REAL pdpsrf(ngrid) ! tendence surf pressure 42 47 REAL pzlev(ngrid,nlay+1) ! altitude at layer boundaries 43 48 REAL pzlay(ngrid,nlay) ! altitude at the middle of the layers 44 49 REAL pt(ngrid,nlay) ! temperature at the middle of the layers (K) 45 REAL pdt(ngrid,nlay) ! tend ance temperature des autres param.50 REAL pdt(ngrid,nlay) ! tendence temperature des autres param. 46 51 47 52 real pq(ngrid,nlay,nq) ! traceur (kg/kg) 48 real pdq(ngrid,nlay,nq) ! tend ance avant condensation (kg/kg.s-1)49 50 REAL tau(ngridmx,naerkind) 51 REAL tauscaling(ngridmx) 52 real rdust(ngridmx,nlay ermx)! Dust geometric mean radius (m)53 real pdq(ngrid,nlay,nq) ! tendence avant condensation (kg/kg.s-1) 54 55 REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point 56 REAL tauscaling(ngridmx) ! Convertion factor for dust amount 57 real rdust(ngridmx,nlay) ! Dust geometric mean radius (m) 53 58 54 59 c Outputs: 55 60 c ------- 56 61 57 real pdqcloud(ngrid,nlay,nq) ! tend ance de la condensation H2O(kg/kg.s-1)62 real pdqcloud(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1) 58 63 real pdqscloud(ngrid,nq) ! flux en surface (kg.m-2.s-1) 59 REAL pdtcloud(ngrid,nlay) ! tend ance temperature due60 ! 64 REAL pdtcloud(ngrid,nlay) ! tendence temperature due 65 ! a la chaleur latente 61 66 62 67 REAL rice(ngrid,nlay) ! Ice mass mean radius (m) … … 64 69 REAL nuice(ngrid,nlay) ! Estimated effective variance 65 70 ! of the size distribution 66 real rsedcloud(ngridmx,nlay ermx) ! Cloud sedimentation radius67 real rhocloud(ngridmx,nlay ermx) ! Cloud density (kg.m-3)71 real rsedcloud(ngridmx,nlay) ! Cloud sedimentation radius 72 real rhocloud(ngridmx,nlay) ! Cloud density (kg.m-3) 68 73 69 74 c local: 70 75 c ------ 71 76 72 INTEGER ig,l 77 ! for sedimentation 78 REAL zq(ngridmx,nlay,nqmx) ! local value of tracers 79 real masse (ngridmx,nlay) ! Layer mass (kg.m-2) 80 real epaisseur (ngridmx,nlay) ! Layer thickness (m) 81 real wq(ngridmx,nlay+1) ! displaced tracer mass (kg.m-2) 82 83 ! for ice radius computation 84 REAL Mo,No 85 REAl ccntyp 86 real beta ! correction for the shape of the ice particles (cf. newsedim) 87 save beta 88 89 ! for time loop 90 INTEGER microstep ! time subsampling step variable 91 INTEGER imicro ! time subsampling for coupled water microphysics & sedimentation 92 SAVE imicro 93 REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation 94 SAVE microtimestep 95 96 ! tendency given by clouds (inside the micro loop) 97 REAL subpdqcloud(ngrid,nlay,nq) ! cf. pdqcloud 98 REAL subpdtcloud(ngrid,nlay) ! cf. pdtcloud 99 100 ! global tendency (clouds+sedim+physics) 101 REAL subpdq(ngrid,nlay,nq) ! cf. pdqcloud 102 REAL subpdt(ngrid,nlay) ! cf. pdtcloud 103 104 REAL CBRT 105 EXTERNAL CBRT 106 107 108 INTEGER iq,ig,l 73 109 LOGICAL,SAVE :: firstcall=.true. 74 110 … … 93 129 write(*,*) "watercloud: igcm_h2o_vap=",igcm_h2o_vap 94 130 write(*,*) " igcm_h2o_ice=",igcm_h2o_ice 131 132 133 write(*,*) "correction for the shape of the ice particles ?" 134 beta=0.75 ! default value 135 call getin("ice_shape",beta) 136 write(*,*) " ice_shape = ",beta 137 138 write(*,*) "time subsampling for microphysic ?" 139 imicro = 1 140 call getin("imicro",imicro) 141 write(*,*)"imicro = ",imicro 142 143 microtimestep = ptimestep/float(imicro) 144 write(*,*)"Physical timestep is",ptimestep 145 write(*,*)"Microphysics timestep is",microtimestep 95 146 96 147 firstcall=.false. 97 148 ENDIF ! of IF (firstcall) 98 99 c Main call to the different cloud schemes: 100 IF (microphys) THEN 101 CALL improvedclouds(ngrid,nlay,ptimestep, 102 & pplev,pplay,pt,pdt, 103 & pq,pdq,pdqcloud,pdqscloud,pdtcloud, 149 150 c-----Initialization 151 subpdq(1:ngrid,1:nlay,1:nq) = 0 152 subpdt(1:ngrid,1:nlay) = 0 153 pdqscloud(1:ngrid,1:nq) = 0 154 zq(1:ngrid,1:nlay,1:nq) = 0 155 156 c-----Computing the different layer properties for clouds sedimentation 157 do l=1,nlay 158 do ig=1, ngrid 159 masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g 160 epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l) 161 enddo 162 enddo 163 164 c------------------------------------------------------------------ 165 c Time subsampling for coupled microphysics and sedimentation 166 c------------------------------------------------------------------ 167 DO microstep=1,imicro 168 169 c------------------------------------------------------------------- 170 c 1. Tendencies: 171 c------------------ 172 173 174 c------ Temperature tendency subpdt 175 ! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt 176 ! If imicro=1 subpdt is the same as pdt 177 DO l=1,nlay 178 DO ig=1,ngrid 179 subpdt(ig,l) = subpdt(ig,l) 180 & + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry 181 ENDDO 182 ENDDO 183 c------ Traceurs tendencies subpdq 184 c------ At each micro timestep we add pdq in order to have a stepped entry 185 IF (microphys) THEN 186 DO l=1,nlay 187 DO ig=1,ngrid 188 subpdq(ig,l,igcm_dust_mass) = 189 & subpdq(ig,l,igcm_dust_mass) 190 & + pdq(ig,l,igcm_dust_mass) 191 subpdq(ig,l,igcm_dust_number) = 192 & subpdq(ig,l,igcm_dust_number) 193 & + pdq(ig,l,igcm_dust_number) 194 subpdq(ig,l,igcm_ccn_mass) = 195 & subpdq(ig,l,igcm_ccn_mass) 196 & + pdq(ig,l,igcm_ccn_mass) 197 subpdq(ig,l,igcm_ccn_number) = 198 & subpdq(ig,l,igcm_ccn_number) 199 & + pdq(ig,l,igcm_ccn_number) 200 ENDDO 201 ENDDO 202 ENDIF 203 DO l=1,nlay 204 DO ig=1,ngrid 205 subpdq(ig,l,igcm_h2o_ice) = 206 & subpdq(ig,l,igcm_h2o_ice) 207 & + pdq(ig,l,igcm_h2o_ice) 208 subpdq(ig,l,igcm_h2o_vap) = 209 & subpdq(ig,l,igcm_h2o_vap) 210 & + pdq(ig,l,igcm_h2o_vap) 211 ENDDO 212 ENDDO 213 214 215 c------------------------------------------------------------------- 216 c 2. Main call to the different cloud schemes: 217 c------------------------------------------------ 218 IF (microphys) THEN 219 CALL improvedclouds(ngrid,nlay,microtimestep, 220 & pplev,pplay,pzlev,pt,subpdt, 221 & pq,subpdq,subpdqcloud,subpdtcloud, 104 222 & nq,tauscaling,rdust,rice,nuice, 105 223 & rsedcloud,rhocloud) 106 ELSE107 CALL simpleclouds(ngrid,nlay,ptimestep,108 & pplev,pplay,pzlev,pzlay,pt, pdt,109 & pq, pdq,pdqcloud,pdqscloud,pdtcloud,224 ELSE 225 CALL simpleclouds(ngrid,nlay,microtimestep, 226 & pplev,pplay,pzlev,pzlay,pt,subpdt, 227 & pq,subpdq,subpdqcloud,subpdtcloud, 110 228 & nq,tau,rice,nuice,rsedcloud) 111 ENDIF 229 ENDIF 230 231 c-------------------------------------------------------------------- 232 c 3. CCN & ice sedimentation: 233 c-------------------------------- 234 ! Sedimentation is done here for water ice and its CCN only 235 ! callsedim in physics is done for dust (and others species if any) 236 237 DO l=1,nlay 238 DO ig=1,ngrid 239 subpdt(ig,l) = 240 & subpdt(ig,l) + subpdtcloud(ig,l) 241 ENDDO 242 ENDDO 243 244 c------- water ice update before sedimentation (radius is done in the cloud scheme routine) 245 DO l=1,nlay 246 DO ig=1, ngrid 247 zq(ig,l,igcm_h2o_ice)= max(1e-30, 248 & pq(ig,l,igcm_h2o_ice) + (subpdq(ig,l,igcm_h2o_ice) 249 & + subpdqcloud(ig,l,igcm_h2o_ice))*microtimestep) 250 ! zq(ig,l,igcm_h2o_vap)= max(1e-30, 251 ! & pq(ig,l,igcm_h2o_vap) + (subpdq(ig,l,igcm_h2o_vap) 252 ! & + subpdqcloud(ig,l,igcm_h2o_vap))*microtimestep) 253 ENDDO 254 ENDDO 255 256 257 c------- CCN update before sedimentation 258 IF (microphys) THEN 259 DO l=1,nlay 260 DO ig=1,ngrid 261 zq(ig,l,igcm_ccn_number)= 262 & pq(ig,l,igcm_ccn_number) + (subpdq(ig,l,igcm_ccn_number) 263 & + subpdqcloud(ig,l,igcm_ccn_number))*microtimestep 264 zq(ig,l,igcm_ccn_number)= max(1e-30, 265 & zq(ig,l,igcm_ccn_number))!*tauscaling(ig)) ! OU pas tauscaling ? 266 zq(ig,l,igcm_ccn_mass)= 267 & pq(ig,l,igcm_ccn_mass) + (subpdq(ig,l,igcm_ccn_mass) 268 & + subpdqcloud(ig,l,igcm_ccn_mass))*microtimestep 269 zq(ig,l,igcm_ccn_mass)=max(1e-30, 270 & zq(ig,l,igcm_ccn_mass))!*tauscaling(ig)) ! OU pas tauscaling ? 271 ENDDO 272 ENDDO 273 ENDIF 274 275 276 277 IF (microphys) THEN 278 279 c------- CCN number sedimentation 280 call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, 281 & microtimestep,pplev,masse,epaisseur,pt,rsedcloud, 282 & rhocloud,zq(1,1,igcm_ccn_number),wq,beta) 283 do ig=1,ngrid 284 ! matters if one would like to know ccn surface deposition 285 pdqscloud(ig,igcm_ccn_number)= 286 & pdqscloud(ig,igcm_ccn_number) + wq(ig,1) 287 enddo 288 289 c------- CCN mass sedimentation 290 call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, 291 & microtimestep,pplev,masse,epaisseur,pt,rsedcloud, 292 & rhocloud,zq(1,1,igcm_ccn_mass),wq,beta) 293 do ig=1,ngrid 294 ! matters if one would like to know ccn surface deposition 295 pdqscloud(ig,igcm_ccn_mass)= 296 & pdqscloud(ig,igcm_ccn_mass) + wq(ig,1) 297 enddo 298 299 c------- Water ice sedimentation -- improved microphys 300 call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, 301 & microtimestep,pplev,masse,epaisseur,pt,rsedcloud, 302 & rhocloud,zq(1,1,igcm_h2o_ice),wq,beta) 303 ELSE 304 305 c------- Water ice sedimentation -- simple microphys 306 call newsedim(ngrid,nlay,ngrid*nlay,1, 307 & microtimestep,pplev,masse,epaisseur,pt,rsedcloud, 308 & rho_q(igcm_h2o_ice),zq(1,1,igcm_h2o_ice),wq,beta) 309 310 ENDIF 311 312 313 c------- Surface ice tendency update 314 DO ig=1,ngrid 315 pdqscloud(ig,igcm_h2o_ice)= 316 & pdqscloud(ig,igcm_h2o_ice) + wq(ig,1) 317 ENDDO 318 319 320 c------------------------------------------------------------------- 321 c 5. Updating tendencies after sedimentation: 322 c----------------------------------------------- 323 324 DO l=1,nlay 325 DO ig=1,ngrid 326 327 subpdq(ig,l,igcm_h2o_ice) = 328 & (zq(ig,l,igcm_h2o_ice) - pq(ig,l,igcm_h2o_ice)) 329 & /microtimestep 330 331 332 subpdq(ig,l,igcm_h2o_vap)=subpdq(ig,l,igcm_h2o_vap) 333 & +subpdqcloud(ig,l,igcm_h2o_vap) 334 335 ENDDO 336 ENDDO 337 338 339 IF (microphys) then 340 DO l=1,nlay 341 DO ig=1,ngrid 342 subpdq(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number) 343 & -pq(ig,l,igcm_ccn_number))/microtimestep 344 subpdq(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass) 345 & -pq(ig,l,igcm_ccn_mass))/microtimestep 346 ENDDO 347 ENDDO 348 ENDIF 349 350 351 352 ENDDO ! of DO microstep=1,imicro 353 354 c------------------------------------------------------------------- 355 c 6. Compute final tendencies after time loop: 356 c------------------------------------------------ 357 358 c------ Whole temperature tendency pdtcloud 359 DO l=1,nlay 360 DO ig=1,ngrid 361 pdtcloud(ig,l) = 362 & subpdt(ig,l)/imicro-pdt(ig,l) 363 ENDDO 364 ENDDO 365 c------ Traceurs tendencies pdqcloud 366 DO iq=1,nq 367 DO l=1,nlay 368 DO ig=1,ngrid 369 pdqcloud(ig,l,iq) = subpdq(ig,l,iq)/imicro 370 & - pdq(ig,l,iq) 371 ENDDO 372 ENDDO 373 ENDDO 374 c------ Traceurs surface tendencies pdqscloud 375 DO iq=1,nq 376 DO ig=1,ngrid 377 pdqscloud(ig,iq) = 378 & pdqscloud(ig,iq)/ptimestep 379 ENDDO 380 ENDDO 381 382 383 384 c------Update the ice particle size "rice" for output or photochemistry. 385 c------It is not used again in the water cycle. 386 IF(scavenging) THEN 387 DO l=1, nlay 388 DO ig=1,ngrid 389 Mo = pq(ig,l,igcm_h2o_ice) 390 & + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep 391 & + (pq(ig,l,igcm_ccn_mass) 392 & + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep) 393 & *tauscaling(ig) + 1.e-30 394 No = (pq(ig,l,igcm_ccn_number) 395 & + pdqcloud(ig,l,igcm_ccn_number)*ptimestep) 396 & *tauscaling(ig) + 1.e-30 397 rhocloud(ig,l) = (pq(ig,l,igcm_h2o_ice) + 398 & pdqcloud(ig,l,igcm_h2o_ice)*ptimestep) 399 & / Mo * rho_ice 400 & + (pq(ig,l,igcm_ccn_mass) 401 & + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep) 402 & *tauscaling(ig)/ Mo * rho_dust 403 rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust) 404 rice(ig,l)=(Mo / No * 0.75 / pi / rhocloud(ig,l))**(1./3.) 405 if ((Mo.lt.1.e-15) .or. (No.le.50)) rice(ig,l) = 1.e-8 406 ENDDO 407 ENDDO 408 ELSE 409 DO l=1,nlay 410 DO ig=1,ngrid 411 ccntyp = 412 & 1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-pzlay(ig,l)/10000.) 413 ccntyp = ccntyp /ccn_factor 414 rice(ig,l)=max( CBRT ( ((pq(ig,l,igcm_h2o_ice) 415 & + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)/rho_ice 416 & +ccntyp*(4./3.)*pi*rdust(ig,l)**3.) 417 & /(ccntyp*4./3.*pi) ), rdust(ig,l)) 418 ENDDO 419 ENDDO 420 ENDIF ! of IF(scavenging) 421 422 !-------------------------------------------------------------- 423 !-------------------------------------------------------------- 424 112 425 113 426 c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005 … … 122 435 end do 123 436 437 c======================================================================= 438 439 !!!!!!!!!! FOR PHOTOCHEMISTRY, REIMPLEMENT output of surfdust/surfice 440 !! if (photochem) then 441 !!c computation of dust and ice surface area (micron2/cm3) 442 !!c for heterogeneous chemistry 443 !! 444 !! do l = 1,nlay 445 !! do ig = 1,ngrid 446 !!c 447 !!c npart: number density of ccn in #/cm3 448 !!c 449 !! npart(ig,l) = 1.e-6*ccn(ig,l) 450 !! $ *masse(ig,l)/epaisseur(ig,l) 451 !!c 452 !!c dust and ice surface area 453 !!c 454 !! surfdust(ig,l) = npart(ig,l)*4.*pi*1.e12*rdust(ig,l)**2 455 !!c 456 !! if (rice(ig,l) .ge. rdust(ig,l)) then 457 !! surfice(ig,l) = npart(ig,l)*4.*pi*1.e12*rice(ig,l)**2 458 !! surfdust(ig,l) = 0. 459 !! else 460 !! surfice(ig,l) = 0. 461 !! end if 462 !! end do ! of do ig=1,ngrid 463 !! end do ! of do l=1,nlay 464 !! end if ! of photochem 465 124 466 RETURN 125 467 END
Note: See TracChangeset
for help on using the changeset viewer.
