Changeset 726
- Timestamp:
- Jul 17, 2012, 5:04:27 PM (12 years ago)
- Location:
- trunk/LMDZ.GENERIC
- Files:
-
- 2 added
- 15 edited
-
README (modified) (1 diff)
-
libf/phystd/aeropacity.F90 (modified) (9 diffs)
-
libf/phystd/aeroptproperties.F90 (modified) (2 diffs)
-
libf/phystd/aerosol_mod.F90 (added)
-
libf/phystd/callcorrk.F90 (modified) (6 diffs)
-
libf/phystd/callkeys.h (modified) (2 diffs)
-
libf/phystd/iniaerosol.F (added)
-
libf/phystd/inifis.F (modified) (4 diffs)
-
libf/phystd/largescale.F (modified) (1 diff)
-
libf/phystd/physdem1.F (modified) (2 diffs)
-
libf/phystd/physiq.F90 (modified) (11 diffs)
-
libf/phystd/radcommon_h.F90 (modified) (3 diffs)
-
libf/phystd/radinc_h.F90 (modified) (1 diff)
-
libf/phystd/rcm1d.F (modified) (3 diffs)
-
libf/phystd/suaer_corrk.F90 (modified) (6 diffs)
-
libf/phystd/tabfi.F (modified) (6 diffs)
-
libf/phystd/turbdiff.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/README
r711 r726 688 688 - Some cleanup in start2archive.F and ini_archive.F to get them to work for 689 689 the generic model (removed some "Martian" specificities). 690 691 == 17/07/2012 == JL for LK 692 - Generalization of aerosol scheme: 693 - any number of aerosols can be used and id numbers are determined consistently by the code. Aerosol order 694 not important anymore. 695 - addition of a module with the id numbers for aerosols (aerosol_mod.F90). 696 - initialization of aerosols id numbers in iniaerosol.F90 697 - compile with -s x where x *must* be equal to the number of aerosols turned on in callphys.def (either by a 698 flag or by dusttau>0 for dust). 699 => may have to erase object files when compiling with s option for the first time. 700 - For no aerosols, run with aeroco2=.true. and aerofixco2=.true (the default distribution for fixed co2 701 aerosols is 1.e-9; can be changed in aeropacity). 702 - If starting from an old start file, recreate start file with the q=0 option in newstart.e. 703 - update callphys.def with aeroXXX and aerofixXXX options (only XXX=co2,h2o supported for 704 now). Dust is activated by setting dusttau>0. See the early mars case in deftank. 705 - To add other aerosols, see Laura Kerber. 706 -
trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90
r716 r726 2 2 aerosol,reffrad,QREFvis3d,QREFir3d,tau_col,cloudfrac,totcloudfrac,clearsky) 3 3 4 use radinc_h, only : naerkind, L_TAUMAX 4 use radinc_h, only : L_TAUMAX,naerkind 5 use aerosol_mod 5 6 6 7 implicit none … … 43 44 #include "comgeomfi.h" 44 45 #include "tracer.h" 46 #include "comvert.h" 45 47 46 48 INTEGER ngrid,nlayer,nq … … 76 78 real topdust, expfactor, zp 77 79 REAL taudusttmp(ngridmx) ! Temporary dust opacity used before scaling 78 80 REAL tauh2so4tmp(ngridmx) ! Temporary h2so4 opacity used before scaling 79 81 ! BENJAMIN MODIFS 80 82 real CLFtot,CLF1,CLF2 … … 100 102 endif 101 103 104 write(*,*) "Tracers found in aeropacity:" 102 105 do iq=1,nqmx 103 106 tracername=noms(iq) 104 107 if (tracername.eq."co2_ice") then 105 108 i_co2ice=iq 109 write(*,*) "i_co2ice=",i_co2ice 110 106 111 endif 107 112 if (tracername.eq."h2o_ice") then 108 113 i_h2oice=iq 114 write(*,*) "i_h2oice=",i_h2oice 109 115 endif 110 116 enddo 111 112 write(*,*) "aeropacity: i_co2ice=",i_co2ice 113 write(*,*) " i_h2oice=",i_h2oice 114 115 if(watercond.and.(.not.aerofixed))then 116 if(naerkind.lt.2)then 117 print*,'Cannot have active H2O clouds with naerkind less than 2!' 118 call abort 119 endif 120 endif 121 122 if(dusttau.gt.0.0)then 123 if(naerkind.lt.3)then 124 print*,'Cannot have active dust with naerkind less than 3!' 125 call abort 126 endif 127 endif 117 118 if (naerkind.ne.0) then 119 print*, "If you would like to use aerosols, make sure any old" 120 print*, "start files are updated in newstart using the option" 121 print*, "q=0" 122 write(*,*) "Aerosols found in aeropacity:" 123 endif 124 125 if (iaero_co2.ne.0) then 126 print*, 'iaero_co2= ',iaero_co2 127 endif 128 if (iaero_h2o.ne.0) then 129 print*,'iaero_h2o= ',iaero_h2o 130 endif 131 if (iaero_dust.ne.0) then 132 print*,'iaero_dust= ',iaero_dust 133 endif 134 if (iaero_h2so4.ne.0) then 135 print*,'iaero_h2so4= ',iaero_h2so4 136 endif 128 137 129 138 firstcall=.false. … … 131 140 132 141 133 DO iaer = 1, naerkind ! Loop on aerosol kind (NOT tracer)134 142 ! --------------------------------------------------------- 135 aerkind: SELECT CASE (iaer) 136 !================================================================== 137 CASE(1) aerkind ! CO2 ice (iaer=1) 138 !================================================================== 139 143 !================================================================== 144 ! CO2 ice aerosols 145 !================================================================== 146 147 if (iaero_co2.ne.0) then 148 iaer=iaero_co2 140 149 ! 1. Initialization 141 150 aerosol(:,:,iaer)=0.0 142 143 151 ! 2. Opacity calculation 144 if (aerofix ed.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed152 if (aerofixco2.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed 145 153 do ig=1, ngrid 146 154 do l=1,nlayer … … 162 170 aerosol(ig,l,iaer) = aerosol0 163 171 ! aerosol(ig,l,iaer) = 0.0 164 172 ! print*, aerosol(ig,l,iaer) 165 173 ! using cloud fraction 166 174 ! aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF)) … … 170 178 ENDDO 171 179 ENDDO 172 end if 173 174 !================================================================== 175 CASE(2) aerkind ! Water ice / liquid (iaer=2) 176 !================================================================== 177 180 end if ! if fixed or varying 181 end if ! if CO2 aerosols 182 !================================================================== 183 ! Water ice / liquid 184 !================================================================== 185 186 if (iaero_h2o.ne.0) then 187 iaer=iaero_h2o 178 188 ! 1. Initialization 179 189 aerosol(:,:,iaer)=0.0 180 181 190 ! 2. Opacity calculation 182 if (aerofix ed.or.(i_h2oice.eq.0).or.clearsky) then191 if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then 183 192 do ig=1,ngrid ! to stop the rad tran bug 184 193 do l=1,nlayer … … 252 261 enddo 253 262 end if 254 255 256 !================================================================== 257 CASE(3) aerkind ! Dust (iaer=3) 258 !================================================================== 263 end if ! End if h2o aerosol 264 265 !================================================================== 266 ! Dust 267 !================================================================== 268 if (iaero_dust.ne.0) then 269 iaer=iaero_dust 270 271 ! 1. Initialization 272 aerosol(:,:,iaer)=0.0 273 274 topdust=30.0 ! km (used to be 10.0 km) LK 275 276 ! 2. Opacity calculation 277 278 ! expfactor=0. 279 DO l=1,nlayer-1 280 DO ig=1,ngrid 281 ! Typical mixing ratio profile 282 283 zp=(preff/pplay(ig,l))**(70./topdust) 284 expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3) 285 286 ! Vertical scaling function 287 aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) & 288 *expfactor 289 290 291 ENDDO 292 ENDDO 293 294 ! Rescaling each layer to reproduce the choosen (or assimilated) 295 ! dust extinction opacity at visible reference wavelength, which 296 ! is scaled to the "preff" reference surface pressure available in comvert.h 297 ! and stored in startfi.nc 298 299 taudusttmp(1:ngrid)=0. 300 DO l=1,nlayer 301 DO ig=1,ngrid 302 taudusttmp(ig) = taudusttmp(ig) & 303 + aerosol(ig,l,iaer) 304 ENDDO 305 ENDDO 306 DO l=1,nlayer-1 307 DO ig=1,ngrid 308 aerosol(ig,l,iaer) = max(1E-20, & 309 dusttau & 310 * pplev(ig,1) / preff & 311 * aerosol(ig,l,iaer) & 312 / taudusttmp(ig)) 313 314 ENDDO 315 ENDDO 316 end if ! If dust aerosol 317 318 !================================================================== 319 ! H2SO4 320 !================================================================== 321 ! added by LK 322 if (iaero_h2so4.ne.0) then 323 iaer=iaero_h2so4 259 324 260 325 ! 1. Initialization 261 326 aerosol(:,:,iaer)=0.0 262 327 263 topdust=10.0 ! km264 265 print*,'WARNING, dust is experimental for Early Mars'266 328 267 329 ! 2. Opacity calculation 268 330 269 do ig=1,ngrid 270 do l=1,nlayer 271 zp=700./pplay(ig,l) 272 aerosol(ig,l,iaer)=(dusttau/700.)*(pplev(ig,l)-pplev(ig,l+1)) & 273 *max( exp(.03*(1.-max(zp,1.))) , 1.E-3 ) & 274 *QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer) 275 ! takes into account particle variation 276 enddo 277 enddo 278 279 !================================================================== 280 END SELECT aerkind 281 ENDDO ! iaer (loop on aerosol kind) 282 331 ! expfactor=0. 332 DO l=1,nlayer-1 333 DO ig=1,ngrid 334 ! Typical mixing ratio profile 335 336 zp=(preff/pplay(ig,l))**(70./30) !emulating topdust 337 expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3) 338 339 ! Vertical scaling function 340 aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) & 341 *expfactor 342 343 344 ENDDO 345 ENDDO 346 tauh2so4tmp(1:ngrid)=0. 347 DO l=1,nlayer 348 DO ig=1,ngrid 349 tauh2so4tmp(ig) = tauh2so4tmp(ig) & 350 + aerosol(ig,l,iaer) 351 ENDDO 352 ENDDO 353 DO l=1,nlayer-1 354 DO ig=1,ngrid 355 aerosol(ig,l,iaer) = max(1E-20, & 356 1 & 357 * pplev(ig,1) / preff & 358 * aerosol(ig,l,iaer) & 359 / tauh2so4tmp(ig)) 360 361 ENDDO 362 ENDDO 363 364 ! 1/700. is assuming a "sulfurtau" of 1 365 ! Sulfur aerosol routine to be improved. 366 ! aerosol0 = & 367 ! ( 0.75 * QREFvis3d(ig,l,iaer) / & 368 ! ( rho_h2so4 * reffrad(ig,l,iaer) ) ) * & 369 ! ( pq(ig,l,i_h2so4) + 1.E-9 ) * & 370 ! ( pplev(ig,l) - pplev(ig,l+1) ) / g 371 ! aerosol0 = max(aerosol0,1.e-9) 372 ! aerosol0 = min(aerosol0,L_TAUMAX) 373 ! aerosol(ig,l,iaer) = aerosol0 374 375 ! ENDDO 376 ! ENDDO 377 end if 378 379 283 380 284 381 ! -------------------------------------------------------------------------- … … 316 413 endif 317 414 end do 318 319 415 return 320 416 end subroutine aeropacity -
trunk/LMDZ.GENERIC/libf/phystd/aeroptproperties.F90
r253 r726 5 5 ! omegaREFvis3d,omegaREFir3d) 6 6 7 use radinc_h, only: L_NSPECTI,L_NSPECTV,n aerkind,nsizemax7 use radinc_h, only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind 8 8 use radcommon_h, only: QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir 9 9 use radcommon_h, only: qrefvis,qrefir,omegarefvis,omegarefir … … 208 208 209 209 if (firstcall) then 210 print*,'Optical prop ertiesof the aerosol are homogenous for:'210 print*,'Optical prop. of the aerosol are homogenous for:' 211 211 print*,'iaer = ',iaer 212 212 endif -
trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
r716 r726 14 14 use ioipsl_getincom 15 15 use gases_h 16 use aerosol_mod 16 17 17 18 implicit none … … 189 190 REAL*8 muvarrad(L_LEVELS) 190 191 191 !======================================================================= 192 193 !=============================================================== 192 194 ! Initialization on first call 193 195 … … 211 213 ! Effective radius and variance of the aerosols 212 214 213 do iaer=1,naerkind 215 214 216 ! these values will change once the microphysics gets to work 215 217 ! UNLESS tracer=.false., in which case we should be working with 216 218 ! a fixed aerosol layer, and be able to define reffrad in a 217 219 ! .def file. To be improved! 218 219 if(iaer.eq.1)then ! CO2 ice 220 do l=1,nlayer 221 do ig=1,ngrid 220 ! Generalization of aerosols added by LK 221 do iaer = 1,naerkind 222 if (iaer.eq.iaero_co2) then ! active CO2 aerosols 223 do l=1,nlayer 224 do ig=1,ngrid 222 225 reffrad(ig,l,iaer) = 1.e-4 223 226 nueffrad(ig,l,iaer) = 0.1 224 enddo 225 enddo 226 endif 227 228 if(iaer.eq.2)then ! H2O ice 229 do l=1,nlayer 230 do ig=1,ngrid 227 enddo 228 enddo 229 endif 230 if (iaer.eq.iaero_h2o) then ! active H2O aerosols 231 do l=1,nlayer 232 do ig=1,ngrid 231 233 reffrad(ig,l,iaer) = 1.e-5 232 234 nueffrad(ig,l,iaer) = 0.1 233 enddo 234 enddo 235 endif 236 237 if(iaer.eq.3)then ! dust 238 do l=1,nlayer 239 do ig=1,ngrid 235 enddo 236 enddo 237 endif 238 if(iaer.eq.iaero_dust)then ! active dust 239 do l=1,nlayer 240 do ig=1,ngrid 240 241 reffrad(ig,l,iaer) = 1.e-5 241 242 nueffrad(ig,l,iaer) = 0.1 242 enddo243 enddo244 endif245 246 if(iaer.gt.3)then247 print*,'Error in callcorrk, naerkind is too high.'248 print*,'The code still needs generalisation to arbitrary'249 print*,'aerosol kinds and number.'250 call abort251 endif252 253 enddo243 enddo 244 enddo 245 endif 246 if(iaer.eq.iaero_h2so4)then ! Active h2so4 aerosols 247 do l=1,nlayer 248 do ig=1,ngrid 249 reffrad(ig,l,iaer) = 1.e-6 250 nueffrad(ig,l,iaer) = 0.1 251 enddo 252 enddo 253 endif 254 enddo 254 255 255 256 print*, "callcorrk: Correlated-k data base folder:",trim(datadir) … … 313 314 314 315 315 316 if(radfixed)then 317 do l=1,nlayer 316 do iaer=1,naerkind 317 318 if(iaer.eq.iaero_co2)then 319 if(radfixed)then 320 do l=1,nlayer 318 321 do ig=1,ngrid 319 reffrad(ig,l, 1) = 5.e-5 ! CO2 ice322 reffrad(ig,l,iaer) = 5.e-5 ! CO2 ice 320 323 enddo 321 enddo 322 print*,'CO2 ice particle size = ',reffrad(1,1,1)/1.e-6,' um' 323 if(naerkind.ge.2)then 324 do l=1,nlayer 325 do ig=1,ngrid 326 !reffrad(ig,l,2) = 2.e-5 ! H2O ice 327 reffrad(ig,l,2) = rad_chaud ! H2O ice 328 329 zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds) 330 zfice = MIN(MAX(zfice,0.0),1.0) 331 reffrad(ig,l,2)= rad_chaud * (1.-zfice) + rad_froid * zfice 332 nueffrad(ig,l,2) = coef_chaud * (1.-zfice) + coef_froid * zfice 333 enddo 334 enddo 335 print*,'H2O ice particle size = ',reffrad(1,1,2)/1.e-6,' um' 336 endif 337 if(naerkind.eq.3)then 338 do l=1,nlayer 339 do ig=1,ngrid 340 reffrad(ig,l,naerkind) = 2.e-6 ! dust 341 enddo 342 enddo 343 print*,'Dust particle size = ',reffrad(1,1,naerkind)/1.e-6,' um' 344 endif 345 if(naerkind.gt.3)then 346 print*,'Code not general enough to deal with naerkind > 3 yet.' 347 call abort 348 endif 349 350 else 351 352 maxrad=0.0 353 !averad=0.0 354 minrad=1.0 355 do l=1,nlayer 356 357 !masse = (pplev(ig,l) - pplev(ig,l+1))/g 324 enddo 325 print*,'CO2 ice particle size =',reffrad(1,1,iaer)/1.e-6,'um' 326 else 327 maxrad=0.0 328 !averad=0.0 329 minrad=1.0 330 do l=1,nlayer 331 332 !masse = (pplev(ig,l) - pplev(ig,l+1))/g 358 333 359 334 do ig=1,ngrid 360 !if(tracer)then361 335 if(tracer.and.igcm_co2_ice.gt.0)then 362 336 363 337 if(igcm_co2_ice.lt.1)then 364 print*,'Tracers but no CO2 ice still seems to be a problem...' 338 print*,'Tracers but no CO2 ice still' 339 print*,'seems to be a problem...' 365 340 print*,'Aborting in callcorrk.' 366 341 stop 367 342 endif 368 369 reffrad(ig,l,1) = CBRT( 3*pq(ig,l,igcm_co2_ice)/ & 343 reffrad(ig,l,iaer) = CBRT( 3*pq(ig,l,igcm_co2_ice)/ & 370 344 (4*Nmix_co2*pi*rho_co2) ) 371 345 endif 372 reffrad(ig,l,1) = max(reffrad(ig,l,1),1.e-6) 373 reffrad(ig,l,1) = min(reffrad(ig,l,1),500.e-6) 374 375 !averad = averad + reffrad(ig,l,1)*area(ig) 376 maxrad = max(reffrad(ig,l,1),maxrad) 377 minrad = min(reffrad(ig,l,1),minrad) 346 reffrad(ig,l,iaer) = max(reffrad(ig,l,iaer),1.e-6) 347 reffrad(ig,l,iaer) = min(reffrad(ig,l,iaer),500.e-6) 348 !averad = averad + reffrad(ig,l,iaer)*area(ig) 349 maxrad = max(reffrad(ig,l,iaer),maxrad) 350 minrad = min(reffrad(ig,l,iaer),minrad) 378 351 enddo 379 enddo 380 if(igcm_co2_ice.gt.0)then 381 print*,'Max. CO2 ice particle size = ',maxrad/1.e-6,' um' 382 print*,'Min. CO2 ice particle size = ',minrad/1.e-6,' um' 352 enddo 353 if(igcm_co2_ice.gt.0)then 354 print*,'Max. CO2 ice particle size = ',maxrad/1.e-6,' um' 355 print*,'Min. CO2 ice particle size = ',minrad/1.e-6,' um' 356 endif 383 357 endif 384 385 if((naerkind.ge.2).and.water)then 386 maxrad=0.0 387 minrad=1.0 358 endif 359 360 361 if(iaer.eq.iaero_h2o)then 362 if(radfixed) then 388 363 do l=1,nlayer 389 364 do ig=1,ngrid 390 reffrad(ig,l,2) = CBRT( 3*pq(ig,l,igcm_h2o_ice)/ &391 (4*Nmix_h2o*pi*rho_ice) )392 reffrad(ig,l,2) = max(reffrad(ig,l,2),1.e-6) 393 reffrad(ig,l,2) = min(reffrad(ig,l,2),100.e-6)394 395 maxrad = max(reffrad(ig,l,2),maxrad)396 minrad = min(reffrad(ig,l,2),minrad)365 !reffrad(ig,l,iaer) = 2.e-5 ! H2O ice 366 reffrad(ig,l,iaer) = rad_chaud ! H2O ice 367 368 zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds) 369 zfice = MIN(MAX(zfice,0.0),1.0) 370 reffrad(ig,l,2)= rad_chaud * (1.-zfice) + rad_froid * zfice 371 nueffrad(ig,l,2) = coef_chaud * (1.-zfice) + coef_froid * zfice 397 372 enddo 398 373 enddo 399 print*,'Max. H2O ice particle size = ',maxrad/1.e-6,' um' 400 print*,'Min. H2O ice particle size = ',minrad/1.e-6,' um' 374 print*,'H2O ice particle size=',reffrad(1,1,iaer)/1.e-6,'um' 375 else 376 if(water)then 377 maxrad=0.0 378 minrad=1.0 379 do l=1,nlayer 380 do ig=1,ngrid 381 reffrad(ig,l,iaer) = CBRT( 3*pq(ig,l,igcm_h2o_ice)/ & 382 (4*Nmix_h2o*pi*rho_ice) ) 383 reffrad(ig,l,iaer) = max(reffrad(ig,l,iaer),1.e-6) 384 reffrad(ig,l,iaer) = min(reffrad(ig,l,iaer),100.e-6) 385 386 maxrad = max(reffrad(ig,l,iaer),maxrad) 387 minrad = min(reffrad(ig,l,iaer),minrad) 388 enddo 389 enddo 390 print*,'Max. H2O ice particle size = ',maxrad/1.e-6,' um' 391 print*,'Min. H2O ice particle size = ',minrad/1.e-6,' um' 392 endif 401 393 endif 402 403 if(naerkind.eq.3)then 394 endif 395 396 if(iaer.eq.iaero_dust)then 404 397 do l=1,nlayer 405 398 do ig=1,ngrid 406 reffrad(ig,l, naerkind) = 2.e-6 ! dust399 reffrad(ig,l,iaer) = 2.e-6 ! dust 407 400 enddo 408 401 enddo 402 print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um' 409 403 endif 410 404 411 endif 405 if(iaer.eq.iaero_h2so4)then 406 do l=1,nlayer 407 do ig=1,ngrid 408 reffrad(ig,l,iaer) = 1.e-6 ! h2so4 409 enddo 410 enddo 411 print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um' 412 endif 413 enddo ! do iaer=1,naerkind 414 ! added by LK 412 415 413 416 … … 523 526 do nw=1,L_NSPECTV 524 527 if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then 525 print*,'Serious problems with qsvaer values in callcorrk' 528 print*,'Serious problems with qsvaer values' 529 print*,'in callcorrk' 526 530 call abort 527 531 endif … … 533 537 do nw=1,L_NSPECTI 534 538 if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then 535 print*,'Serious problems with qsiaer values in callcorrk' 539 print*,'Serious problems with qsiaer values' 540 print*,'in callcorrk' 536 541 call abort 537 542 endif -
trunk/LMDZ.GENERIC/libf/phystd/callkeys.h
r716 r726 34 34 & , sedimentation,water,watercond,waterrain & 35 35 & , rainthreshold & 36 & , aerofixed & 36 & , aeroco2,aeroh2o,aeroh2so4 & 37 & , aerofixco2,aerofixh2o & 37 38 & , hydrology & 38 39 & , sourceevol & … … 72 73 logical sedimentation 73 74 logical water,watercond,waterrain 74 logical aerofixed 75 logical aeroco2,aeroh2o,aeroh2so4 76 logical aerofixco2,aerofixh2o 75 77 logical hydrology 76 78 logical sourceevol -
trunk/LMDZ.GENERIC/libf/phystd/inifis.F
r716 r726 367 367 ! TRACERS: 368 368 369 write(*,*)"Fixed / inactive aerosol distribution?"370 aerofixed=.true. ! default value371 call getin("aerofixed",aerofixed)372 write(*,*)" aerofixed = ",aerofixed373 374 369 write(*,*)"Varying H2O cloud fraction?" 375 370 CLFvarying=.false. ! default value … … 392 387 write(*,*)" Nmix_h2o = ",Nmix_h2o 393 388 389 ! write(*,*)"Number of radiatively active aerosols:" 390 ! naerkind=0. ! default value 391 ! call getin("naerkind",naerkind) 392 ! write(*,*)" naerkind = ",naerkind 393 394 394 write(*,*)"Opacity of dust (if used):" 395 395 dusttau=0. ! default value … … 397 397 write(*,*)" dusttau = ",dusttau 398 398 399 ! Test of incompatibility: 400 ! if less than 3 aerosols, then dusttau should = 0 401 if ((naerkind.lt.3).and.dusttau.gt.0.) then 402 print*,'Need naer>2 for radiatively active dust!' 403 stop 404 endif 399 write(*,*)"Radiatively active CO2 aerosols?" 400 aeroco2=.false. ! default value 401 call getin("aeroco2",aeroco2) 402 write(*,*)" aeroco2 = ",aeroco2 403 404 write(*,*)"Fixed CO2 aerosol distribution?" 405 aerofixco2=.false. ! default value 406 call getin("aerofixco2",aerofixco2) 407 write(*,*)" aerofixco2 = ",aerofixco2 408 409 write(*,*)"Radiatively active water ice?" 410 aeroh2o=.false. ! default value 411 call getin("aeroh2o",aeroh2o) 412 write(*,*)" aeroh2o = ",aeroh2o 413 414 write(*,*)"Fixed H2O aerosol distribution?" 415 aerofixh2o=.false. ! default value 416 call getin("aerofixh2o",aerofixh2o) 417 write(*,*)" aerofixh2o = ",aerofixh2o 418 419 write(*,*)"Radiatively active sulfuric acid aersols?" 420 aeroh2so4=.false. ! default value 421 call getin("aeroh2so4",aeroh2so4) 422 write(*,*)" aeroh2so4 = ",aeroh2so4 405 423 406 424 write(*,*)"Cloud pressure level (with kastprof only):" … … 426 444 427 445 ! Test of incompatibility: 428 ! if no tracers, then aerofixed should be true 429 if ((.not.tracer).and.(.not.aerofixed)) then 430 print*,'if tracers are off, aerofixed must be ON!' 446 ! if no tracers, then the variable aerosols should be turned off 447 if ((.not.tracer).and.(.not.aerofixco2)) then 448 print*,'if tracers are off, you must turn aerofixco2 on.' 449 print*, '(i.e., you must fix the aerosol distribution)' 450 print*, 'This option was formerly known as aerofixed' 431 451 stop 432 452 endif 453 if ((.not.tracer).and.(.not.aerofixh2o)) then 454 print*,'if tracers are off, you must turn aerofixh2o on.' 455 print*, '(i.e., you must fix the aerosol distribution)' 456 print*, 'This option was formerly known as aerofixed' 457 stop 458 endif 459 433 460 434 461 ! Test of incompatibility: -
trunk/LMDZ.GENERIC/libf/phystd/largescale.F
r253 r726 111 111 112 112 ! for varying particle size in rad tran and (possibly) sedimentation 113 if(.not.aerofixed)then 113 if(aeroh2o.and.(.not.aerofixh2o)) then 114 114 115 reffH2O(i,k) = CBRT( 3*zcond(i)/( 4*Nmix_h2o*pi*rho_ice)) 115 116 reffH2O(i,k) = max(reffH2O(i,k),1.e-8) -
trunk/LMDZ.GENERIC/libf/phystd/physdem1.F
r590 r726 5 5 . cloudfrac,totcloudfrac,hice) 6 6 7 use radcommon_h, only: tauvis8 7 9 8 implicit none … … 215 214 tab_cntrl(33) = dtemisice(1) ! time scale for snow metamorphism (north) 216 215 tab_cntrl(34) = dtemisice(2) ! time scale for snow metamorphism (south) 217 218 c dust aerosol properties219 tab_cntrl(27) = tauvis ! mean visible optical depth220 216 221 217 tab_cntrl(28) = 0. -
trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
r716 r726 7 7 pdu,pdv,pdt,pdq,pdpsrf,tracerdyn) 8 8 9 use radinc_h, only : naerkind,L_NSPECTI,L_NSPECTV9 use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind 10 10 use watercommon_h, only : RLVTT 11 11 use gases_h 12 12 use radcommon_h, only: sigma 13 use aerosol_mod 13 14 implicit none 14 15 … … 389 390 ! included by RW for variable H2O particle sizes 390 391 real reffH2O(ngridmx,nlayermx) 391 real reffcol(ngridmx, 2)392 real reffcol(ngridmx,naerkind) 392 393 393 394 ! included by RW for sourceevol … … 420 421 zdtsw(:,:) = 0.0 421 422 zdtlw(:,:) = 0.0 423 424 ! initialize aerosol indexes 425 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 426 call iniaerosol() 427 422 428 423 429 ! initialize tracer names, indexes and properties … … 427 433 call initracer() 428 434 endif ! end tracer 435 436 ! 429 437 430 438 ! read startfi (initial parameters) … … 811 819 dEtotLW = cpp*SUM(massarea(:,:)*zdtlw(:,:))/totarea 812 820 dEtotsSW = SUM(fluxsurf_sw(:)*(1.-albedo(:))*area(:))/totarea 813 dEtotsLW = SUM( fluxsurf_lw(:)*emis(:)*area(:)-zplanck(:))/totarea821 dEtotsLW = SUM((fluxsurf_lw(:)*emis(:)-zplanck(:))*area(:))/totarea 814 822 dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:) 815 823 dEzRadlw(:,:)=cpp*mass(:,:)*zdtlw(:,:) … … 1045 1053 call abort 1046 1054 endif 1047 1055 print*, 'we are in co2cond!!!!!' 1048 1056 call condense_cloud(ngrid,nlayer,nq,ptimestep, & 1049 1057 capcal,pplay,pplev,tsurf,pt, & … … 1161 1169 pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+dqcldlscale(ig,l) 1162 1170 pdt(ig,l) = pdt(ig,l)+dtlscale(ig,l) 1163 1164 if( .not.aerofixed)then1165 reffrad(ig,l, 2)=reffH2O(ig,l)1171 1172 if(aeroh2o.and.(.not.aerofixh2o)) then 1173 reffrad(ig,l,iaero_h2o)=reffH2O(ig,l) 1166 1174 endif 1167 1175 … … 1541 1549 enddo 1542 1550 close(38) 1543 print*,'As testradtimes enabled, exiting physics on first call' 1551 print*,'As testradtimes enabled,' 1552 print*,'exiting physics on first call' 1544 1553 call abort 1545 1554 endif … … 1559 1568 enddo 1560 1569 1561 ! not generalised for arbitrary aerosols yet!!!1570 ! Generalised for arbitrary aerosols now. (LK) 1562 1571 reffcol(:,:)=0.0 1563 1572 do ig=1,ngrid 1564 1573 do l=1,nlayer 1565 if(co2cond )then1566 reffcol(ig, 1) = reffcol(ig,1) + zq(ig,l,igcm_co2_ice) * &1567 reffrad(ig,l, 1) * &1574 if(co2cond.and.(iaero_co2.ne.0)) then 1575 reffcol(ig,iaero_co2) = reffcol(ig,iaero_co2) + zq(ig,l,igcm_co2_ice) * & 1576 reffrad(ig,l,iaero_co2) * & 1568 1577 (pplev(ig,l) - pplev(ig,l+1)) / g 1569 1578 endif 1570 if(water )then1571 reffcol(ig, 2) = reffcol(ig,2) + zq(ig,l,igcm_h2o_ice) * &1572 reffrad(ig,l, 2) * &1579 if(water.and.(iaero_h2o.ne.0)) then 1580 reffcol(ig,iaero_h2o) = reffcol(ig,iaero_h2o) + zq(ig,l,igcm_h2o_ice) * & 1581 reffrad(ig,l,iaero_h2o) * & 1573 1582 (pplev(ig,l) - pplev(ig,l+1)) / g 1574 1583 endif … … 1732 1741 call wstats(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col', & 1733 1742 'kg m^-2',2,qcol(1,iq) ) 1734 call wstats(ngridmx,trim(noms(iq))//'_reff', &1735 trim(noms(iq))//'_reff', &1736 'm',3,reffrad(1,1,iq))1743 ! call wstats(ngridmx,trim(noms(iq))//'_reff', & 1744 ! trim(noms(iq))//'_reff', & 1745 ! 'm',3,reffrad(1,1,iq)) 1737 1746 end do 1738 1747 if (water) then … … 1809 1818 1810 1819 ! Output aerosols 1811 if (igcm_co2_ice.ne.0) call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,1)) 1812 if (igcm_h2o_ice.ne.0) call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(1,1,2)) 1813 if (igcm_co2_ice.ne.0) call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,1)) 1814 if (igcm_h2o_ice.ne.0) call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,2)) 1820 if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) & 1821 call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3, & 1822 reffrad(1,1,iaero_co2)) 1823 if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) & 1824 call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3, & 1825 reffrad(1,1,iaero_h2o)) 1826 if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) & 1827 call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol', & 1828 'um kg m^-2',2,reffcol(1,iaero_co2)) 1829 if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) & 1830 call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol', & 1831 'um kg m^-2',2,reffcol(1,iaero_h2o)) 1815 1832 1816 1833 ! Output tracers -
trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90
r484 r726 46 46 ! ~~~~~~~~~ 47 47 ! 48 ! tauvis: dust optical depth at reference wavelength ("longrefvis" set49 ! in dimradmars.h : typically longrefvis = 0.67E-6 m, as measured by Viking )50 !51 48 ! For the "naerkind" kind of aerosol radiative properties : 52 49 ! QVISsQREF : Qext / Qext("longrefvis") … … 61 58 ! omegaIR : mean single scattering albedo 62 59 ! gIR : mean assymetry factor 63 64 60 65 61 REAL*8 BWNI(L_NSPECTI+1), WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI) … … 78 74 real*8 tgasmin, tgasmax 79 75 80 real tauvis81 76 real QVISsQREF(L_NSPECTV,naerkind,nsizemax) 82 77 real omegavis(L_NSPECTV,naerkind,nsizemax) -
trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90
r716 r726 5 5 #include "dimensions.h" 6 6 #include "bands.h" 7 #include "scatterers.h"7 #include "scatterers.h" 8 8 9 9 !====================================================================== -
trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F
r650 r726 1 1 program rcm1d 2 2 3 use radcommon_h, only: tauvis4 3 ! to use 'getin' 5 4 use ioipsl_getincom … … 469 468 write(*,*)" specOLR = ",specOLR 470 469 471 !!! AS12: what shall we do with this?472 c Proprietes des poussiere aerosol473 c --------------------------------474 tauvis=0.2 ! default value for tauvis475 print *,'Reference dust opacity at 700 Pa ?'476 call getin("tauvis",tauvis)477 write(*,*) " tauvis = ",tauvis478 479 470 c 480 471 c pour le schema d'ondes de gravite … … 709 700 ! write(103,*) 'Ls=',zls*180./pi 710 701 ! write(103,*) 'Lat=', lati(1)*180./pi 711 ! write(103,*) 'Tau=', tauvis/700*psurf712 702 ! write(103,*) 'RunEnd - Atmos. Temp. File' 713 703 ! write(103,*) 'RunEnd - Atmos. Temp. File' 714 704 ! write(104,*) 'Ls=',zls*180./pi 715 705 ! write(104,*) 'Lat=', lati(1) 716 ! write(104,*) 'Tau=', tauvis/700*psurf717 706 ! write(104,*) 'RunEnd - Atmos. Temp. File' 718 707 ENDIF -
trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90
r374 r726 2 2 3 3 ! inputs 4 use radinc_h, only: L_NSPECTI,L_NSPECTV,n aerkind,nsizemax,iim,jjm4 use radinc_h, only: L_NSPECTI,L_NSPECTV,nsizemax,iim,jjm,naerkind 5 5 use radcommon_h, only: blamv,blami,lamrefir,lamrefvis 6 6 use datafile_mod, only: datadir … … 10 10 use radcommon_h, only: radiustab,nsize,tstellar 11 11 use radcommon_h, only: qrefvis,qrefir,omegarefvis,omegarefir 12 use aerosol_mod 12 13 13 14 implicit none … … 44 45 45 46 #include "callkeys.h" 46 47 47 ! Optical properties (read in external ASCII files) 48 48 INTEGER :: nwvl ! Number of wavelengths in … … 101 101 ! Please also choose the reference wavelengths of each aerosol 102 102 103 forgetit=.true. 104 105 ! naerkind=1 103 do iaer=1,naerkind 104 if (iaer.eq.iaero_co2)then 105 forgetit=.true. 106 print*, 'naerkind= co2', iaer 107 106 108 ! visible 107 if(forgetit)then108 file_id(1,1) = 'optprop_co2_vis_ff.dat' ! Francois' values109 else110 file_id(1,1) = 'optprop_co2ice_vis_n50.dat'111 endif112 lamrefvis(1)=1.5E-6 ! 1.5um OMEGA/MEx ???109 if(forgetit)then 110 file_id(iaer,1) = 'optprop_co2_vis_ff.dat' ! Francois' values 111 else 112 file_id(iaer,1) = 'optprop_co2ice_vis_n50.dat' 113 endif 114 lamrefvis(iaer)=1.5E-6 ! 1.5um OMEGA/MEx ??? 113 115 114 116 ! infrared 115 if(forgetit)then116 file_id(1,2) = 'optprop_co2_ir_ff.dat' ! Francois' values117 else118 file_id(1,2) = 'optprop_co2ice_ir_n50.dat'119 endif120 lamrefir(1)=12.1E-6 ! 825cm-1 TES/MGS ???121 117 if(forgetit)then 118 file_id(iaer,2) = 'optprop_co2_ir_ff.dat' ! Francois' values 119 else 120 file_id(iaer,2) = 'optprop_co2ice_ir_n50.dat' 121 endif 122 lamrefir(iaer)=12.1E-6 ! 825cm-1 TES/MGS ??? 123 endif ! CO2 aerosols 122 124 ! NOTE: these lamref's are currently for dust, not CO2 ice. 123 125 ! JB thinks this shouldn't matter too much, but it needs to be 124 126 ! fixed / decided for the final version. 125 127 126 IF (naerkind .gt. 1) THEN 127 ! naerkind=2 128 if (iaer.eq.iaero_h2o) then 129 print*, 'naerkind= h2o', iaer 130 128 131 ! visible 129 file_id( 2,1) = 'optprop_icevis_n50.dat'130 lamrefvis( 2)=1.5E-6 ! 1.5um OMEGA/MEx132 file_id(iaer,1) = 'optprop_icevis_n50.dat' 133 lamrefvis(iaer)=1.5E-6 ! 1.5um OMEGA/MEx 131 134 ! infrared 132 file_id(2,2) = 'optprop_iceir_n50.dat' 133 lamrefir(2)=12.1E-6 ! 825cm-1 TES/MGS 134 ENDIF 135 136 IF (naerkind .eq. 3) THEN 137 ! naerkind=3 135 file_id(iaer,2) = 'optprop_iceir_n50.dat' 136 lamrefir(iaer)=12.1E-6 ! 825cm-1 TES/MGS 137 endif 138 139 if (iaer.eq.iaero_dust) then 140 print*, 'naerkind= dust', iaer 141 138 142 ! visible 139 file_id( naerkind,1) = 'optprop_dustvis_n50.dat'140 !lamrefvis( 3)=1.5E-6 ! 1.5um OMEGA/MEx141 lamrefvis( naerkind)=0.67e-6143 file_id(iaer,1) = 'optprop_dustvis_n50.dat' 144 !lamrefvis(iaer)=1.5E-6 ! 1.5um OMEGA/MEx 145 lamrefvis(iaer)=0.67e-6 142 146 ! infrared 143 file_id(naerkind,2) = 'optprop_dustir_n50.dat' 144 !lamrefir(3)=12.1E-6 ! 825cm-1 TES/MGS 145 lamrefir(naerkind)=9.3E-6 146 147 ENDIF 148 149 IF (naerkind .gt. 3) THEN 150 print*,'naerkind = ',naerkind 151 print*,'but we only have data for 3 types, exiting.' 152 call abort 153 ENDIF 147 file_id(iaer,2) = 'optprop_dustir_n50.dat' 148 !lamrefir(iaer)=12.1E-6 ! 825cm-1 TES/MGS 149 lamrefir(iaer)=9.3E-6 150 endif 151 152 if (iaer.eq.iaero_h2so4) then 153 print*, 'naerkind= h2so4', iaer 154 155 ! visible 156 file_id(iaer,1) = 'optprop_h2so4vis_n20.dat' 157 !lamrefir(iaer)= doesn't exist? 158 lamrefvis(iaer)=1.5E-6 ! no idea, must find 159 ! infrared 160 file_id(iaer,2) = 'optprop_h2so4ir_n20.dat' 161 !lamrefir(iaer)=12.1E-6 ! 825cm-1 TES/MGS 162 lamrefir(iaer)=9.3E-6 ! no idea, must find 163 ! added by LK 164 endif 165 166 enddo 154 167 155 168 !------------------------------------------------------------------ … … 282 295 283 296 jfile = jfile+1 284 IF ((idomain.NE. 1).OR.(iaer.NE.1)) THEN297 IF ((idomain.NE.iaero_co2).OR.(iaer.NE.iaero_co2)) THEN 285 298 endwhile = .true. 286 299 ENDIF … … 300 313 301 314 ! 1.5 If Francois' data, convert wvl to metres 302 if(forgetit)then 315 if(iaer.eq.iaero_co2)then 316 if(forgetit)then 303 317 ! print*,'please sort out forgetit for naerkind>1' 304 if(iaer.eq.1)then305 318 do iwvl=1,nwvl 306 319 wvl(iwvl)=wvl(iwvl)*1.e-6 -
trunk/LMDZ.GENERIC/libf/phystd/tabfi.F
r649 r726 45 45 ! to use 'getin' 46 46 use ioipsl_getincom 47 use radcommon_h, only: tauvis48 47 49 48 implicit none … … 138 137 emisice(2) = tab_cntrl(tab0+25) 139 138 emissiv = tab_cntrl(tab0+26) 140 tauvis = tab_cntrl(tab0+27) ! dust opt. depth vis.141 139 iceradius(1)= tab_cntrl(tab0+31) ! mean scat radius of CO2 snow (north) 142 140 iceradius(2)= tab_cntrl(tab0+32) ! mean scat radius of CO2 snow (south) … … 200 198 write(*,5) '(33) dtemisice(1)',tab_cntrl(tab0+33),dtemisice(1) 201 199 write(*,5) '(34) dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2) 202 203 write(*,5) '(27) tauvis',tab_cntrl(tab0+27),tauvis204 200 205 201 write(*,5) '(35) volcapa',tab_cntrl(tab0+35),volcapa … … 228 224 write(*,*) '(31 et 32) iceradius : mean scat radius of CO2 snow' 229 225 write(*,*) '(33 et 34) dtemisice : time scale for snow', 230 & ' metamorphism' 231 write(*,*) '(27) tauvis : mean dust vis. reference ', 232 & ' opacity' 233 write(*,*) '(35) volcapa : soil volumetric heat capacity' 226 & 'metamorphism' 227 write(*,*) '(35) volcapa : soil volumetric heat capacity' 234 228 write(*,*) '(18) obliquit : planet obliquity (deg)' 235 229 write(*,*) '(17) peri_day : periastron date (sols since Ls=0)' … … 365 359 write(*,*) ' dtemisice(2) (new value):',dtemisice(2) 366 360 367 else if (modif(1:len_trim(modif)) .eq. 'tauvis') then368 write(*,*) 'current value:',tauvis369 write(*,*) 'enter new value:'370 114 read(*,*,iostat=ierr) tauvis371 if(ierr.ne.0) goto 114372 write(*,*)373 write(*,*) ' tauvis (new value):',tauvis374 375 361 else if (modif(1:len_trim(modif)) .eq. 'obliquit') then 376 362 write(*,*) 'current value:',obliquit … … 535 521 write(*,5) '(34) dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2) 536 522 537 write(*,5) '(27) tauvis',tab_cntrl(tab0+27),tauvis538 539 523 write(*,5) '(35) volcapa',tab_cntrl(tab0+35),volcapa 540 524 -
trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90
r650 r726 137 137 138 138 if(water)then 139 ivap=igcm_h2o_vap 140 iliq=igcm_h2o_ice 141 iliq_surf=igcm_h2o_vap 142 iice_surf=igcm_h2o_ice ! simply to make the code legible 143 ! to be generalised later 139 ivap=igcm_h2o_vap 140 iliq=igcm_h2o_ice 141 iliq_surf=igcm_h2o_vap 142 iice_surf=igcm_h2o_ice ! simply to make the code legible 143 ! to be generalised 144 else 145 ivap=0 146 iliq=0 147 iliq_surf=0 148 iice_surf=0 ! simply to make the code legible 144 149 endif 145 150 sensibFlux(:)=0.
Note: See TracChangeset
for help on using the changeset viewer.
