Changeset 2622
- Timestamp:
- Feb 16, 2022, 4:26:05 PM (3 years ago)
- Location:
- trunk/LMDZ.VENUS/libf/phyvenus
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.VENUS/libf/phyvenus/aeropacity.F90
r2560 r2622 2 2 aerosol,reffrad,nueffrad,QREFvis3d,tau_col) 3 3 4 use radinc_h, only : L_TAUMAX,naerkind4 use radinc_h, only : naerkind 5 5 use aerosol_mod 6 6 7 7 implicit none 8 8 #include "YOMCST.h" … … 49 49 REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth 50 50 51 real aerosol0, obs_tau_col_aurora, pm 52 real pcloud_deck, cloud_slope 53 54 real dp_strato(ngrid) 55 real dp_tropo(ngrid) 56 real dp_layer(ngrid) 57 58 INTEGER l,ig,iq,iaer,ia,ll 51 INTEGER l,ig,iaer,ll 59 52 60 53 LOGICAL,SAVE :: firstcall=.true. 61 54 !$OMP THREADPRIVATE(firstcall) 62 REAL CBRT 63 EXTERNAL CBRT 64 65 ! for fixed dust profiles 66 real topdust, expfactor, zp 67 REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling 68 REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling 69 70 real CLFtot 55 character*80 abort_message 56 character*20 modname 71 57 72 58 ! for venus clouds 73 real :: p_bot2,p_bot,p_top,p_top2,h_bot2,h_bot,h_top,h_top2 74 real :: mode_dens,h_lay,nmode1,nmode2,nmode2p,nmode3,nmodeuv 75 59 real,ALLOCATABLE,SAVE :: pbot2(:,:),pbot(:,:),ptop(:,:),ptop2(:,:) 60 real,ALLOCATABLE,SAVE :: hbot2(:,:),hbot(:,:),htop(:,:),htop2(:,:) 61 real,ALLOCATABLE,SAVE :: nmode(:,:) 62 !$OMP THREADPRIVATE(pbot2,pbot,ptop,ptop2) 63 !$OMP THREADPRIVATE(hbot2,hbot,htop,htop2) 64 !$OMP THREADPRIVATE(nmode) 65 real :: mode_dens,hlay 66 67 ! --------------------------------------------------------- 68 ! --------------------------------------------------------- 76 69 ! First call only 77 70 IF (firstcall) THEN 78 79 ! identify tracers 80 write(*,*) "Active aerosols found in aeropacity:" 81 82 if (iaero_venus1.ne.0) then 71 72 modname = 'aeropacity' 73 74 ! verify tracers 75 write(*,*) "Active aerosols found in aeropacity, designed for Venus" 76 77 if (iaero_venus1.eq.1) then 83 78 print*,'iaero_venus1= ',iaero_venus1 84 endif 85 if (iaero_venus2.ne.0) then 79 else 80 abort_message='iaero_venus1 is not 1...' 81 call abort_physic(modname,abort_message,1) 82 endif 83 if (iaero_venus2.eq.2) then 86 84 print*,'iaero_venus2= ',iaero_venus2 87 endif 88 if (iaero_venus2p.ne.0) then 85 else 86 abort_message='iaero_venus2 is not 2...' 87 call abort_physic(modname,abort_message,1) 88 endif 89 if (iaero_venus2p.eq.3) then 89 90 print*,'iaero_venus2p= ',iaero_venus2p 90 endif 91 if (iaero_venus3.ne.0) then 91 else 92 abort_message='iaero_venus2p is not 3...' 93 call abort_physic(modname,abort_message,1) 94 endif 95 if (iaero_venus3.eq.4) then 92 96 print*,'iaero_venus3= ',iaero_venus3 93 endif 94 if (iaero_venusUV.ne.0) then 97 else 98 abort_message='iaero_venus3 is not 4...' 99 call abort_physic(modname,abort_message,1) 100 endif 101 if (iaero_venusUV.eq.5) then 95 102 print*,'iaero_venusUV= ',iaero_venusUV 96 endif 103 else 104 abort_message='iaero_venus1UV is not 5...' 105 call abort_physic(modname,abort_message,1) 106 endif 107 108 ! cloud model 109 allocate(pbot2(ngrid,naerkind),pbot(ngrid,naerkind)) 110 allocate(ptop(ngrid,naerkind),ptop2(ngrid,naerkind)) 111 allocate(hbot2(ngrid,naerkind),hbot(ngrid,naerkind)) 112 allocate(htop(ngrid,naerkind),htop2(ngrid,naerkind)) 113 allocate(nmode(ngrid,naerkind)) 114 call cloud_haus16(ngrid,pbot2,pbot,hbot2,hbot,ptop,ptop2,htop,htop2,nmode) 97 115 98 116 firstcall=.false. 99 117 ENDIF ! of IF (firstcall) 100 101 102 ! --------------------------------------------------------- 103 104 !================================================================== 105 ! Venus clouds (4 modes) 106 ! S. Lebonnois (jan 2016) 107 !================================================================== 108 ! distributions from Haus et al, 2016 109 ! mode 1 2 2p 3 110 ! r (microns) 0.30 1.05 1.40 3.65 111 ! sigma 1.56 1.29 1.23 1.28 112 ! reff (microns) 0.49 1.23 1.56 4.25 113 ! nueff 0.21 0.067 0.044 0.063 114 ! (nueff=exp(ln^2 sigma)-1) 115 ! 116 ! p_bot <=> zb ; p_top <=> zb+zc ; h_bot <=> Hlo ; h_top <=> Hup 117 ! p<p_top: N(i+1)=N(i)*(p(i+1)/p(i))**(h_lay(i)/h_top) 118 ! h_lay=R(T(i)+T(i+1))/2/g (in m) 119 ! R=8.31/mu (mu in kg/mol) 120 ! p>p_bot: N(i-1)=N(i)*(p(i)/p(i-1))**(h_lay(i)/h_bot) 121 ! N is in m-3 122 ! 123 ! dTau = Qext*[pi*reff**2*exp(-3*ln(1+nueff))]*N*h_lay*(-dp)/p 124 125 ! Mode 1 126 if (iaero_venus1 .ne.0) then 127 iaer=iaero_venus1 128 129 ! 1. Initialization 130 aerosol(1:ngrid,1:nlayer,iaer)=0.0 131 ! two scales heights for below and above max density (cf Haus et al 16) 132 p_bot2 = 2.0e5 ! Pa 2.0e5 133 h_bot2 = 5.0e3 ! m 134 p_bot = 1.2e5 ! 1.2e5 135 h_bot = 1.0e3 136 nmode1 = 1.935e8 ! m-3 137 p_top = 1.0e4 138 h_top = 3.5e3 139 p_top2 = 4.5e2 140 h_top2 = 2.0e3 141 142 ! 2. Opacity calculation 143 144 DO ig=1,ngrid 118 ! --------------------------------------------------------- 119 ! --------------------------------------------------------- 120 121 aerosol(:,:,:) = 0. 122 tau_col(:) = 0. 123 124 do ig=1,ngrid 125 do iaer=1,naerkind 126 127 ! Opacity calculation 128 145 129 ! determine the approximate middle of the mode layer 146 130 ll=1 147 131 DO l=1,nlayer-1 ! high p to low p 148 IF (pplay(ig,l) .gt. (p _top+p_bot)/2) ll=l132 IF (pplay(ig,l) .gt. (ptop(ig,iaer)+pbot(ig,iaer))/2) ll=l 149 133 ENDDO 150 ! from there go down and upfor profile N151 mode_dens = nmode 1! m-3134 ! from there go DOWN for profile N 135 mode_dens = nmode(ig,iaer) ! m-3 152 136 DO l=ll,1,-1 153 h _lay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG154 IF (pplay(ig,l) .le. p _bot) THEN155 mode_dens = nmode 1! m-3156 ELSEIF (pplay(ig,l) .gt. p _bot .and. pplay(ig,l) .le. p_bot2) THEN157 mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h _lay/h_bot)158 ELSEIF (pplay(ig,l) .gt. p _bot2) THEN159 mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h _lay/h_bot2)137 hlay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG 138 IF (pplay(ig,l) .le. pbot(ig,iaer)) THEN 139 mode_dens = nmode(ig,iaer) ! m-3 140 ELSEIF (pplay(ig,l) .gt. pbot(ig,iaer) .and. pplay(ig,l) .le. pbot2(ig,iaer)) THEN 141 mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(hlay/hbot(ig,iaer)) 142 ELSEIF (pplay(ig,l) .gt. pbot2(ig,iaer)) THEN 143 mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(hlay/hbot2(ig,iaer)) 160 144 ENDIF 161 145 mode_dens=max(1.e-30,mode_dens) 162 aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)* & 163 RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 164 mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l) 165 ! aerosol(ig,l,iaer) = mode_dens 146 if (iaer.lt.5) then 147 aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)* & 148 RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 149 mode_dens*hlay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l) 150 else ! for UV absorber 151 ! normalized to 0.35 microns (peak of absorption) 152 aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens 153 endif 166 154 ENDDO 167 mode_dens = nmode1 ! m-3 155 ! ... then UP for profile N 156 mode_dens = nmode(ig,iaer) ! m-3 168 157 DO l=ll+1,nlayer-1 169 h _lay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG170 IF (pplay(ig,l) .gt. p _top) THEN171 mode_dens = nmode 1! m-3172 ELSEIF (pplay(ig,l) .le. p _top .and. pplay(ig,l) .gt. p_top2) THEN173 mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h _lay/h_top)174 ELSEIF (pplay(ig,l) .le. p _top2) THEN175 mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h _lay/h_top2)158 hlay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG 159 IF (pplay(ig,l) .gt. ptop(ig,iaer)) THEN 160 mode_dens = nmode(ig,iaer) ! m-3 161 ELSEIF (pplay(ig,l) .le. ptop(ig,iaer) .and. pplay(ig,l) .gt. ptop2(ig,iaer)) THEN 162 mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(hlay/htop(ig,iaer)) 163 ELSEIF (pplay(ig,l) .le. ptop2(ig,iaer)) THEN 164 mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(hlay/htop2(ig,iaer)) 176 165 ENDIF 177 166 mode_dens=max(1.e-30,mode_dens) 178 aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)* & 179 RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 180 mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l) 181 ! aerosol(ig,l,iaer) = mode_dens 167 if (iaer.lt.5) then 168 aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)* & 169 RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 170 mode_dens*hlay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l) 171 else ! for UV absorber 172 ! normalized to 0.35 microns (peak of absorption) 173 aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens 174 endif 182 175 ENDDO 183 ENDDO 184 185 end if ! mode 1 186 187 ! Mode 2 188 if (iaero_venus2 .ne.0) then 189 iaer=iaero_venus2 190 191 ! 1. Initialization 192 aerosol(1:ngrid,1:nlayer,iaer)=0.0 193 p_bot = 1.0e4 194 h_bot = 3.0e3 195 nmode2 = 1.00e8 ! m-3 196 p_top = 8.0e3 197 h_top = 3.5e3 198 p_top2 = 4.5e2 199 h_top2 = 2.0e3 200 201 ! 2. Opacity calculation 202 203 DO ig=1,ngrid 204 ! determine the approximate middle of the mode layer 205 ll=1 206 DO l=1,nlayer-1 ! high p to low p 207 IF (pplay(ig,l) .gt. (p_top+p_bot)/2) ll=l 208 ENDDO 209 ! from there go down and up for profile N 210 mode_dens = nmode2 ! m-3 211 DO l=ll,1,-1 212 h_lay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG 213 IF (pplay(ig,l) .le. p_bot) THEN 214 mode_dens = nmode2 ! m-3 215 ELSEIF (pplay(ig,l) .gt. p_bot) THEN 216 mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h_lay/h_bot) 217 ENDIF 218 mode_dens=max(1.e-30,mode_dens) 219 aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)* & 220 RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 221 mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l) 222 ! aerosol(ig,l,iaer) = mode_dens 223 ENDDO 224 mode_dens = nmode2 ! m-3 225 DO l=ll+1,nlayer-1 226 h_lay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG 227 IF (pplay(ig,l) .gt. p_top) THEN 228 mode_dens = nmode2 ! m-3 229 ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. p_top2) THEN 230 mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top) 231 ELSEIF (pplay(ig,l) .le. p_top2) THEN 232 mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top2) 233 ENDIF 234 mode_dens=max(1.e-30,mode_dens) 235 aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)* & 236 RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 237 mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l) 238 ! aerosol(ig,l,iaer) = mode_dens 239 ENDDO 240 ENDDO 241 242 end if ! mode 2 243 244 ! Mode 2p 245 if (iaero_venus2p .ne.0) then 246 iaer=iaero_venus2p 247 248 ! 1. Initialization 249 aerosol(1:ngrid,1:nlayer,iaer)=0.0 250 p_bot = 1.2e5 ! 1.2e5 251 h_bot = 0.1e3 252 nmode2p = 5.0e7 ! m-3 253 p_top = 2.3e4 254 h_top = 1.0e3 255 256 ! 2. Opacity calculation 257 258 DO ig=1,ngrid 259 ! determine the approximate middle of the mode layer 260 ll=1 261 DO l=1,nlayer-1 ! high p to low p 262 IF (pplay(ig,l) .gt. (p_top+p_bot)/2) ll=l 263 ENDDO 264 ! from there go down and up for profile N 265 mode_dens = nmode2p ! m-3 266 DO l=ll,1,-1 267 h_lay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG 268 IF (pplay(ig,l) .le. p_bot) THEN 269 mode_dens = nmode2p ! m-3 270 ELSEIF (pplay(ig,l) .gt. p_bot) THEN 271 mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h_lay/h_bot) 272 ENDIF 273 mode_dens=max(1.e-30,mode_dens) 274 aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)* & 275 RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 276 mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l) 277 ! aerosol(ig,l,iaer) = mode_dens 278 ENDDO 279 mode_dens = nmode2p ! m-3 280 DO l=ll+1,nlayer-1 281 h_lay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG 282 IF (pplay(ig,l) .gt. p_top) THEN 283 mode_dens = nmode2p ! m-3 284 ELSEIF (pplay(ig,l) .le. p_top) THEN 285 mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top) 286 ENDIF 287 mode_dens=max(1.e-30,mode_dens) 288 aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)* & 289 RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 290 mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l) 291 ! aerosol(ig,l,iaer) = mode_dens 292 ENDDO 293 ENDDO 294 295 end if ! mode 2p 296 297 ! Mode 3 298 if (iaero_venus3 .ne.0) then 299 iaer=iaero_venus3 300 301 ! 1. Initialization 302 aerosol(1:ngrid,1:nlayer,iaer)=0.0 303 p_bot = 1.2e5 ! 1.2e5 304 h_bot = 0.5e3 305 nmode3 = 1.4e7 ! m-3 306 p_top = 3.9e4 307 h_top = 1.0e3 308 309 ! 2. Opacity calculation 310 311 DO ig=1,ngrid 312 ! determine the approximate middle of the mode layer 313 ll=1 314 DO l=1,nlayer-1 ! high p to low p 315 IF (pplay(ig,l) .gt. (p_top+p_bot)/2) ll=l 316 ENDDO 317 ! from there go down and up for profile N 318 mode_dens = nmode3 ! m-3 319 DO l=ll,1,-1 320 h_lay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG 321 IF (pplay(ig,l) .le. p_bot) THEN 322 mode_dens = nmode3 ! m-3 323 ELSEIF (pplay(ig,l) .gt. p_bot) THEN 324 mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h_lay/h_bot) 325 ENDIF 326 mode_dens=max(1.e-30,mode_dens) 327 aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)* & 328 RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 329 mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l) 330 ! aerosol(ig,l,iaer) = mode_dens 331 ENDDO 332 mode_dens = nmode3 ! m-3 333 DO l=ll+1,nlayer-1 334 h_lay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG 335 IF (pplay(ig,l) .gt. p_top) THEN 336 mode_dens = nmode3 ! m-3 337 ELSEIF (pplay(ig,l) .le. p_top) THEN 338 mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top) 339 ENDIF 340 mode_dens=max(1.e-30,mode_dens) 341 aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)* & 342 RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 343 mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l) 344 ! aerosol(ig,l,iaer) = mode_dens 345 ENDDO 346 ENDDO 347 348 end if ! mode 3 349 350 ! UV absorber 351 if (iaero_venusUV .ne.0) then 352 iaer=iaero_venusUV 353 354 ! 1. Initialization 355 aerosol(1:ngrid,1:nlayer,iaer)=0.0 356 p_bot = 3.3e4 ! 58 km 357 h_bot = 1.0e3 358 nmodeuv = 1.00e7 !m-3 359 p_top = 3.7e3 ! 70 km 360 h_top = 1.0e3 361 362 ! 2. Opacity calculation 363 364 DO ig=1,ngrid 365 ! determine the approximate middle of the mode layer 366 ll=1 367 DO l=1,nlayer-1 ! high p to low p 368 IF (pplay(ig,l) .gt. (p_top+p_bot)/2) ll=l 369 ENDDO 370 ! from there go down and up for profile N 371 mode_dens = nmodeuv ! m-3 372 DO l=ll,1,-1 373 h_lay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG 374 IF (pplay(ig,l) .le. p_bot) THEN 375 mode_dens = nmodeuv ! m-3 376 ELSEIF (pplay(ig,l) .gt. p_bot) THEN 377 mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h_lay/h_bot) 378 ENDIF 379 mode_dens=max(1.e-30,mode_dens) 380 ! normalized to 0.35 microns (peak of absorption) 381 aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens 382 ENDDO 383 mode_dens = nmodeuv ! m-3 384 DO l=ll+1,nlayer-1 385 h_lay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG 386 IF (pplay(ig,l) .gt. p_top) THEN 387 mode_dens = nmodeuv ! m-3 388 ELSEIF (pplay(ig,l) .le. p_top) THEN 389 mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top) 390 ENDIF 391 mode_dens=max(1.e-30,mode_dens) 392 ! normalized to 0.35 microns (peak of absorption) 393 aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens 394 ENDDO 395 ENDDO 396 397 ! 3. Re-normalize to Haus et al 2015 total column optical depth 398 tau_col(:)=0.0 399 DO l=1,nlayer 400 DO ig=1,ngrid 176 177 if (iaer.eq.5) then ! for UV absorber 178 DO l=1,nlayer 401 179 tau_col(ig) = tau_col(ig) & 402 180 + aerosol(ig,l,iaer) 403 181 ENDDO 404 ENDDO 405 DO ig=1,ngrid 406 DO l=1,nlayer-1 182 DO l=1,nlayer 407 183 aerosol(ig,l,iaer)=aerosol(ig,l,iaer)*0.205/tau_col(ig) 408 ENDDO 409 ENDDO 410 411 end if ! UV absorber 184 ENDDO 185 endif 186 187 enddo 188 enddo 412 189 413 190 !================================================================== … … 429 206 430 207 tau_col(:)=0.0 431 do i aer = 1, naerkind208 do ig=1,ngrid 432 209 do l=1,nlayer 433 do i g=1,ngrid210 do iaer = 1, naerkind 434 211 tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer) 435 212 end do … … 459 236 endif 460 237 end do 238 461 239 return 462 240 end subroutine aeropacity 241 242 ! -------------------------------------------------------------------------- 243 ! -------------------------------------------------------------------------- 244 subroutine cloud_haus16(ngrid,pbot2,pbot,hbot2,hbot,ptop,ptop2,htop,htop2,nmode) 245 246 !================================================================== 247 ! Venus clouds (4 modes) 248 ! S. Lebonnois (jan 2016) 249 !================================================================== 250 ! distributions from Haus et al, 2016 251 ! ------------------------------------ 252 ! 253 ! mode 1 2 2p 3 254 ! r (microns) 0.30 1.05 1.40 3.65 255 ! sigma 1.56 1.29 1.23 1.28 256 ! reff (microns) 0.49 1.23 1.56 4.25 257 ! nueff 0.21 0.067 0.044 0.063 258 ! (nueff=exp(ln^2 sigma)-1) 259 ! 260 ! pbot <=> zb ; ptop <=> zb+zc ; hbot <=> Hlo ; htop <=> Hup 261 ! p<ptop: N(i+1)=N(i)*(p(i+1)/p(i))**(hlay(i)/htop) 262 ! hlay=R(T(i)+T(i+1))/2/g (in m) 263 ! R=8.31/mu (mu in kg/mol) 264 ! p>pbot: N(i-1)=N(i)*(p(i)/p(i-1))**(hlay(i)/hbot) 265 ! N is in m-3 266 ! 267 ! values in each mode below from Table 1 and text (p.186) 268 ! 269 ! dTau = Qext*[pi*reff**2*exp(-3*ln(1+nueff))]*N*hlay*(-dp)/p 270 ! 271 ! Latitudinal dependence (values from tables 3 and 4): 272 ! mode factor MF12(lat) for modes 1 and 2 273 ! mode 2' unchanged 274 ! mode factor MF3(lat) for mode 3 275 ! + for mode 2 only : pbot(lat), ptop(lat), htop(lat), htop2(lat) 276 !================================================================== 277 278 use radinc_h, only : naerkind 279 USE geometry_mod, ONLY: latitude_deg 280 implicit none 281 282 INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns 283 real,INTENT(out) :: pbot2(ngrid,naerkind),pbot(ngrid,naerkind) 284 real,INTENT(out) :: ptop(ngrid,naerkind),ptop2(ngrid,naerkind) 285 real,INTENT(out) :: hbot2(ngrid,naerkind),hbot(ngrid,naerkind) 286 real,INTENT(out) :: htop(ngrid,naerkind),htop2(ngrid,naerkind) 287 real,INTENT(out) :: nmode(ngrid,naerkind) 463 288 289 ! For latitudinal dependence 290 integer :: i,j,ilat 291 real :: latit,factlat 292 integer,parameter :: nlat_H16 = 16 293 real :: lat_H16(nlat_H16) 294 real :: MF12(nlat_H16),MF3(nlat_H16) 295 real :: pbotM2(nlat_H16),ptopM2(nlat_H16),htopM2(nlat_H16) 296 297 !------------------------- 298 !! Equatorial model 299 !------------------------- 300 ! initialization 301 pbot2(:,:) = 1.e8 302 pbot(:,:) = 1.e8 303 ptop2(:,:) = 0. 304 ptop(:,:) = 0. 305 hbot2(:,:) = 1. 306 hbot(:,:) = 1. 307 htop2(:,:) = 1. 308 htop(:,:) = 1. 309 nmode(:,:) = 0. 310 311 ! table of values(lat,mode) 312 313 ! Mode 1 314 pbot2(:,1) = 2.0e5 ! Pa 315 hbot2(:,1) = 5.0e3 ! m 316 pbot(:,1) = 1.2e5 317 hbot(:,1) = 1.0e3 318 nmode(:,1) = 1.935e8 ! m-3 319 ptop(:,1) = 1.0e4 320 htop(:,1) = 3.5e3 321 ptop2(:,1) = 4.5e2 322 htop2(:,1) = 2.0e3 323 324 ! Mode 2 325 pbot(:,2) = 1.0e4 326 hbot(:,2) = 3.0e3 327 nmode(:,2) = 1.00e8 ! m-3 328 ptop(:,2) = 8.0e3 329 htop(:,2) = 3.5e3 330 ptop2(:,2) = 4.5e2 331 htop2(:,2) = 2.0e3 332 333 ! Mode 2p 334 pbot(:,3) = 1.2e5 335 hbot(:,3) = 0.1e3 336 nmode(:,3) = 5.0e7 ! m-3 337 ptop(:,3) = 2.3e4 338 htop(:,3) = 1.0e3 339 340 ! Mode 3 341 pbot(:,4) = 1.2e5 342 hbot(:,4) = 0.5e3 343 nmode(:,4) = 1.4e7 ! m-3 344 ptop(:,4) = 3.9e4 345 htop(:,4) = 1.0e3 346 347 ! UV absorber 348 pbot(:,5) = 3.3e4 ! 58 km 349 hbot(:,5) = 1.0e3 350 nmode(:,5) = 1.00e7 !m-3 351 ptop(:,5) = 3.7e3 ! 70 km 352 htop(:,5) = 1.0e3 353 354 !---------------------------- 355 !! Latitudinal variations 356 !---------------------------- 357 lat_H16(:) = (/0.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,90./) 358 MF12(:) = (/0.98,0.98,0.99,1.00,0.98,0.94,0.86,0.81,0.73,0.67,0.64,0.61,0.59,0.47,0.36,0.36/) 359 MF3(:) = (/1.30,1.30,1.26,1.23,1.17,1.13,1.06,1.03,1.04,1.09,1.22,1.51,1.82,2.02,2.09,2.09/) 360 pbotM2(:) = (/1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,9.4e3, & 361 9.2e3,9.1e3,9.6e3,1.14e4,1.21e4,1.25e4,1.25e4/) 362 ptopM2(:) = (/8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,7.7e3, & 363 7.5e3,7.5e3,8.0e3,9.2e3,1.0e4,1.06e4,1.06e4/) 364 htopM2(:) = (/3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.4e3, & 365 3.2e3,2.6e3,2.0e3,2.0e3,0.6e3,0.5e3,0.5e3/) 366 367 do i=1,ngrid 368 latit = abs(latitude_deg(i)) 369 ilat=2 370 do j=2,nlat_H16 371 if (latit.gt.lat_H16(j-1)) ilat=j 372 enddo 373 factlat = (lat_H16(ilat)-latit)/(lat_H16(ilat)-lat_H16(ilat-1)) 374 pbot(i,2) = pbotM2(ilat)*(1.-factlat) + pbotM2(ilat-1)*factlat 375 ptop(i,2) = ptopM2(ilat)*(1.-factlat) + ptopM2(ilat-1)*factlat 376 htop(i,2) = htopM2(ilat)*(1.-factlat) + htopM2(ilat-1)*factlat 377 htop2(i,2) = min(htop2(i,2),htop(i,2)) 378 nmode(i,1) = nmode(i,1)*(MF12(ilat)*(1.-factlat)+MF12(ilat-1)*factlat) 379 nmode(i,2) = nmode(i,2)*(MF12(ilat)*(1.-factlat)+MF12(ilat-1)*factlat) 380 nmode(i,3) = nmode(i,3)*(MF3(ilat) *(1.-factlat)+MF3(ilat-1) *factlat) 381 enddo 382 383 return 384 end subroutine cloud_haus16 -
trunk/LMDZ.VENUS/libf/phyvenus/diffusion.h
r2198 r2622 11 11 parameter (Pdiff=1.e-1) ! pressure below which diffusion is computed 12 12 parameter (tdiffmin=5.d0) 13 parameter (dzres=0. 1d0) ! grid resolution (km) for diffusion13 parameter (dzres=0.2d0) ! grid resolution (km) for diffusion 14 14 -
trunk/LMDZ.VENUS/libf/phyvenus/hrtherm.F
r2464 r2622 76 76 xabsi(3,i) = rm(i,i_o) 77 77 xabsi(8,i) = rm(i,i_n2) 78 xabsi(11,i) 78 xabsi(11,i) = rm(i,i_co) 79 79 80 80 c xabsi(6,i) = rm(i,i_h2o2) -
trunk/LMDZ.VENUS/libf/phyvenus/nonoro_gwd_ran_mod.F90
r2580 r2622 96 96 ! generic : kmin=1/sqrt(DxDy) Dx and Dy horizontal grid 97 97 98 ! REAL, parameter:: cmax = 61. ! Max horizontal absolute phase velocity TestGW6 98 !---------------------------------------- 99 ! VCD 1.1 tuning 100 ! REAL, parameter:: cmax = 61. ! Max horizontal absolute phase velocity 99 101 !---------------------------------------- 100 102 ! VCD 2.0 tuning 101 REAL, parameter:: cmax = 301. ! Max horizontal absolute phase velocity TestGW6102 !---------------------------------------- 103 103 REAL, parameter:: cmax = 61. ! Max horizontal absolute phase velocity 104 !---------------------------------------- 105 104 106 ! generic : cmax=zonal wind at launch 105 107 REAL, parameter:: cmin = 1. ! Min horizontal absolute phase velocity … … 126 128 127 129 ! 0.3.2 Parameters of waves dissipations 128 REAL, parameter:: sat = 0.85 ! saturation parameter 129 REAL, parameter:: rdiss = 0.1 ! coefficient of dissipation 130 !---------------------------------------- 131 ! VCD 1.1 tuning 132 ! REAL, parameter:: sat = 0.85 ! saturation parameter 133 ! REAL, parameter:: rdiss = 0.1 ! coefficient of dissipation 134 !---------------------------------------- 135 ! VCD 2.0 tuning 136 REAL, parameter:: sat = 0.6 ! saturation parameter 137 REAL, parameter:: rdiss = 8.e-4 ! coefficient of dissipation 138 !---------------------------------------- 130 139 REAL, parameter:: zoisec = 1.e-8 ! security for intrinsic freguency 131 140 132 141 ! 0.3.3 Background flow at 1/2 levels and vertical coordinate 133 142 REAL H0bis(ngrid, nlayer) ! characteristic height of the atmosphere 143 REAL min_k(ngrid) ! min(kstar,kmin) 134 144 REAL, save:: H0 ! characteristic height of the atmosphere 135 REAL, parameter:: pr = 5e5 ! Reference pressure [Pa] ! VENUS: cloud layer 145 REAL, save:: MDR ! characteristic mass density 146 !---------------------------------------- 147 ! VCD 1.1 tuning 148 ! REAL, parameter:: pr = 5e5 ! Reference pressure [Pa] ! VENUS: cloud layer 149 !---------------------------------------- 150 ! VCD 2.0 tuning 151 REAL, parameter:: pr = 5e4 ! Reference pressure [Pa] ! VENUS: cloud layer 152 !---------------------------------------- 136 153 REAL, parameter:: tr = 300. ! Reference temperature [K] ! VENUS: cloud layer 137 154 REAL zh(ngrid, nlayer + 1) ! Log-pressure altitude (constant H0) … … 156 173 ! Characteristic vertical scale height 157 174 H0 = RD * tr / RG 175 ! Characteristic mass density at launch altitude 176 MDR = pr / ( RD * tr ) 158 177 ! Control 159 178 if (deltat .LT. dtime) THEN … … 214 233 uh(:, nlayer + 1) = uu(:, nlayer) 215 234 vh(:, nlayer + 1) = vv(:, nlayer) 235 236 ! TN+GG April/June 2020 237 ! "Individual waves are not supposed 238 ! to occupy the entire domain, but 239 ! only a faction of it" Lott+2012 240 ! minimum value of k between kmin and kstar 241 DO ii = 1, ngrid 242 kstar = RPI / sqrt(cell_area(ii)) 243 min_k(ii) = max(kmin,kstar) 244 end DO 216 245 217 246 ! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY … … 235 264 ! horizontal wavenumber amplitude 236 265 ! zk(jw, ii) = kmin + (kmax - kmin) * ran_num_2 237 ! TN+GG April/June 2020 238 ! "Individual waves are not supposed 239 ! to occupy the entire domain, but 240 ! only a faction of it" Lott+2012 241 kstar = RPI / sqrt(cell_area(ii)) 242 zk(jw, ii) = max(kmin,kstar) + (kmax - max(kmin,kstar)) * ran_num_2 266 zk(jw, ii) = min_k(ii) + (kmax - min_k(ii)) * ran_num_2 243 267 244 268 ! horizontal phase speed … … 320 344 ! Intrinsic frequency 321 345 intr_freq_p(jw, :) = zo(jw, :) - zk(jw, :) * cos(zp(jw, :)) * uh(:, ll + 1) & 322 - zk(jw, :) * sin(zp(jw, :)) * vh(:, ll + 1)346 - zk(jw, :) * sin(zp(jw, :)) * vh(:, ll + 1) 323 347 ! Vertical velocity in flux formulation 324 348 wwp(jw, :) = min( & … … 336 360 * abs(intr_freq_p(jw, :))**3 / bv(:, ll + 1) & 337 361 * exp(-zh(:, ll + 1) / H0) & 338 * sat**2 * kmin**2 / zk(jw, :)**4) 362 !! Correction for VCD 2.0 363 ! * sat**2 * kmin**2 / zk(jw, :)**4) 364 * sat**2 * min_k(:)**2 * MDR / zk(jw, :)**4) 339 365 end DO 340 366 -
trunk/LMDZ.VENUS/libf/phyvenus/photochemistry_venus.F90
r2580 r2622 1569 1569 !-- TuneC 1570 1570 ! VCD 1.1 tuning 1571 1572 1571 ! v_phot(65:78,4) = v_phot(65:78,4)*10. 1572 ! v_phot(52:59,3) = v_phot(52:59,3)*5. 1573 1573 !-- 1574 1574 !-- TuneE 1575 1575 ! VCD 2.0 tuning 1576 !v_phot(65:90,4) = v_phot(65:90,4)*10.1576 v_phot(65:90,4) = v_phot(65:90,4)*10. 1577 1577 !-- 1578 1578 ! v_phot(:,4) = v_phot(:,4)*10. -
trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F
r2580 r2622 1449 1449 mmean(:,:) = 0. 1450 1450 do iq = 1,nqmax - nmicro 1451 mmean(:,:) = mmean(:,:)+tr_seri(:,:,iq)*m_tr(iq) 1452 rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K 1451 mmean(:,:) = mmean(:,:)+tr_seri(:,:,iq)/m_tr(iq) 1453 1452 enddo 1453 mmean(:,:) = 1./mmean(:,:) 1454 rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K 1454 1455 endif 1455 1456 … … 1812 1813 ! Obsolete 1813 1814 ! but used for VCD 1.1 1814 call flott_gwd_ran(klon,klev,dtime,pplay,zn2, 1815 e t_seri, u_seri, v_seri, paprs(klon/2+1,:), 1815 ! call flott_gwd_ran(klon,klev,dtime,pplay,zn2, 1816 ! e t_seri, u_seri, v_seri, paprs(klon/2+1,:), 1817 ! o zustrhi,zvstrhi, 1818 ! o d_t_hin, d_u_hin, d_v_hin) 1819 1820 ! New routine based on Generic 1821 ! used after VCD 1.1, for VCD 2.0 1822 call nonoro_gwd_ran(klon,klev,dtime,pplay,zn2,presnivs, 1823 e t_seri, u_seri, v_seri, 1816 1824 o zustrhi,zvstrhi, 1817 1825 o d_t_hin, d_u_hin, d_v_hin) 1818 1819 ! New routine based on Generic1820 ! used after VCD 1.11821 ! call nonoro_gwd_ran(klon,klev,dtime,pplay,zn2,presnivs,1822 ! e t_seri, u_seri, v_seri,1823 ! o zustrhi,zvstrhi,1824 ! o d_t_hin, d_u_hin, d_v_hin)1825 1826 1826 1827 c ajout des tendances -
trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F
r2580 r2622 74 74 if (debutphy) then 75 75 76 ! !--- Adjustment of Helium amount77 ! 78 ! trac(:,:,i_he)=trac(:,:,i_he)*20.79 ! 80 ! !---76 !--- Adjustment of Helium amount 77 ! if (i_he/=0) then 78 ! trac(:,:,i_he)=trac(:,:,i_he)/2. 79 ! endif 80 !--- 81 81 82 82 !------------------------------------------------------------------- … … 86 86 87 87 !!! in this reinitialisation, trac is VOLUME mixing ratio 88 print*, "Tracers are re-initialised" 89 trac(:,:,:) = 1.0e-30 90 91 if ((i_ocs /= 0) .and. (i_co /= 0) .and. (i_hcl /= 0) 92 $ .and. (i_so2 /= 0) .and. (i_h2o /= 0) .and. (i_n2/ = 0) 93 $ .and. (i_co2 /= 0)) then 94 95 trac(:,1:22,i_ocs) = 3.e-6 96 trac(:,1:22,i_co) = 25.e-6 97 trac(:,:,i_hcl) = 0.4e-6 98 trac(:,1:22,i_so2) = 7.e-6 99 trac(:,1:22,i_h2o) = 30.e-6 100 trac(:,:,i_n2) = 0.35e-1 88 ! ONLY SO2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 89 ! convert mass to volume mixing ratio 90 do iq = 1,nqmax - nmicro 91 trac(:,:,iq) = trac(:,:,iq)*mmean(:,:)/m_tr(iq) 92 end do 93 print*, "SO2 is re-initialised" 94 if (i_so2 /= 0) then 95 trac(:,1:22,i_so2) = 10.e-6 96 97 ! AL TRACERS!!!!!!!!!!!!!!!!!!!!!!!!!!!! 98 ! print*, "Tracers are re-initialised" 99 ! trac(:,:,:) = 1.0e-30 100 101 ! if ((i_ocs /= 0) .and. (i_co /= 0) .and. (i_hcl /= 0) 102 ! $ .and. (i_so2 /= 0) .and. (i_h2o /= 0) .and. (i_n2/ = 0) 103 ! $ .and. (i_co2 /= 0)) then 104 105 ! trac(:,1:22,i_ocs) = 3.e-6 106 ! trac(:,1:22,i_co) = 25.e-6 107 ! trac(:,:,i_hcl) = 0.4e-6 108 ! trac(:,1:22,i_so2) = 7.e-6 109 ! trac(:,1:22,i_h2o) = 30.e-6 110 ! trac(:,:,i_n2) = 0.35e-1 111 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 101 112 102 113 ! ensure that sum of mixing ratios = 1 … … 124 135 do iq = 1,nqmax - nmicro 125 136 mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq) 126 rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K127 137 enddo 138 rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K 128 139 129 140 ! convert volume to mass mixing ratio … … 293 304 mmean(:,:) = 0. 294 305 do iq = 1,nqmax - nmicro 295 mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq) 296 rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K 306 mmean(:,:) = mmean(:,:)+ztrac(:,:,iq)*m_tr(iq) 297 307 enddo 298 308 rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K 309 310 !=================================================================== 311 ! convert volume to mass mixing ratio / then tendencies in mmr 312 !=================================================================== 299 313 ! gas phase 300 314
Note: See TracChangeset
for help on using the changeset viewer.