Changeset 508
- Timestamp:
- Jan 26, 2012, 3:09:06 PM (13 years ago)
- Location:
- trunk/LMDZ.MARS/libf/phymars
- Files:
-
- 2 edited
-
calltherm_interface.F90 (modified) (11 diffs)
-
thermcell_main_mars.F90 (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90
r499 r508 15 15 #include "dimphys.h" 16 16 #include "comcstfi.h" 17 #include "tracer.h" 17 18 18 19 !-------------------------------------------------------- … … 61 62 REAL ztla(ngridmx,nlayermx) 62 63 REAL q_therm(ngridmx,nlayermx), pq_therm(ngridmx,nlayermx,nqmx) 63 REAL dq_therm(ngridmx,nlayermx), dq_thermdown(ngridmx,nlayermx)64 64 REAL q2_therm(ngridmx,nlayermx), dq2_therm(ngridmx,nlayermx) 65 65 REAL lmax_real(ngridmx) … … 100 100 REAL pbl_teta(ngridmx),dteta(ngridmx,nlayermx) 101 101 102 !--------------------------------------------------------- 103 102 !-------------------------------------------------------- 103 ! Theta_m 104 !-------------------------------------------------------- 105 106 INTEGER ico2 107 SAVE ico2 104 108 105 109 ! ********************************************************************** … … 115 119 q2_therm(:,:)=0. 116 120 dq2_therm(:,:)=0. 117 dq_therm(:,:)=0.118 dq_thermdown(:,:)=0.119 121 ztla(:,:)=0. 120 122 pbl_dtke(:,:)=0. … … 165 167 166 168 ! ********************************************************************** 169 ! Polar night mixing : theta_m 170 ! ********************************************************************** 171 172 if(firstcall) then 173 ico2=0 174 if (tracer) then 175 ! Prepare Special treatment if one of the tracers is CO2 gas 176 do iq=1,nqmx 177 if (noms(iq).eq."co2") then 178 ico2=iq 179 end if 180 enddo 181 endif 182 endif !of if firstcall 183 184 185 ! ********************************************************************** 167 186 ! ********************************************************************** 168 187 ! ********************************************************************** … … 204 223 zheatFlux(:,:)=0. 205 224 zheatFlux_down(:,:)=0. 206 !zbuoyancyOut(:,:)=0.207 !zbuoyancyEst(:,:)=0.225 zbuoyancyOut(:,:)=0. 226 zbuoyancyEst(:,:)=0. 208 227 zzw2(:,:)=0. 209 228 zmax(:)=0. … … 233 252 ! d_v_the(:,:)=d_v_the(:,:)*fact 234 253 ! dq2_the(:,:)=dq2_the(:,:)*fact 235 ! if (nqmx.ne. 0) then236 ! d_q_the(:,:,:)=d_q_the(:,:,:)*fact237 !endif254 if (ico2 .ne. 0) then 255 d_q_the(:,:,ico2)=d_q_the(:,:,ico2)*fact 256 endif 238 257 239 258 zmaxth(:)=zmaxth(:)+zmax(:)*fact … … 250 269 heatFlux_down(:,:)=heatFlux_down(:,:) & 251 270 & +zheatFlux_down(:,:)*fact 252 !buoyancyOut(:,:)=buoyancyOut(:,:) &253 !& +zbuoyancyOut(:,:)*fact254 !buoyancyEst(:,:)=buoyancyEst(:,:) &255 !& +zbuoyancyEst(:,:)*fact271 buoyancyOut(:,:)=buoyancyOut(:,:) & 272 & +zbuoyancyOut(:,:)*fact 273 buoyancyEst(:,:)=buoyancyEst(:,:) & 274 & +zbuoyancyEst(:,:)*fact 256 275 257 276 zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact … … 262 281 ! d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:) 263 282 ! d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:) 264 ! d_q_ajs(:,:,:)=d_q_ajs(:,:,:)+d_q_the(:,:,:) 283 if (ico2 .ne. 0) then 284 d_q_ajs(:,:,ico2)=d_q_ajs(:,:,ico2)+d_q_the(:,:,ico2) 285 endif 265 286 ! dq2_therm(:,:)=dq2_therm(:,:)+dq2_the(:,:) 266 287 ! incrementation des variables meteo … … 269 290 ! zu(:,:) = zu(:,:) + d_u_the(:,:) 270 291 ! zv(:,:) = zv(:,:) + d_v_the(:,:) 271 ! pq_therm(:,:,:) = pq_therm(:,:,:) + d_q_the(:,:,:) 292 if (ico2 .ne. 0) then 293 pq_therm(:,:,ico2) = & 294 & pq_therm(:,:,ico2) + d_q_the(:,:,ico2)*ptimestep 295 endif 272 296 ! q2_therm(:,:) = q2_therm(:,:) + dq2_therm(:,:) 273 297 … … 303 327 modname='tracer' 304 328 DO iq=1,nqmx 329 if (iq .ne. ico2) then 305 330 call thermcell_dqup(ngridmx,nlayermx,ptimestep & 306 331 & ,fm_therm,entr_therm,detr_therm, & 307 332 & masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq),modname,zdz) 308 333 endif 309 334 ENDDO 310 335 endif -
trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90
r496 r508 20 20 #include "dimphys.h" 21 21 #include "comcstfi.h" 22 #include "tracer.h" 23 #include "callkeys.h" 22 24 23 25 ! ============== INPUTS ============== … … 110 112 REAL wmaxa(ngridmx) 111 113 REAL zdz,zbuoy(ngridmx,nlayermx),zw2m 112 LOGICAL active (ngridmx),activetmp(ngridmx)114 LOGICAL activecell(ngridmx),activetmp(ngridmx) 113 115 REAL a1,b1,ae,be,ad,bd,fdfu 114 116 INTEGER tic … … 163 165 ! ========================================= 164 166 167 ! ============== Theta_M Variables ======== 168 169 INTEGER ico2 170 SAVE ico2 171 REAL m_co2, m_noco2, A , B 172 SAVE A, B 173 LOGICAL firstcall 174 save firstcall 175 data firstcall/.true./ 176 REAL zhc(ngridmx,nlayermx) 177 REAL zdzfull(ngridmx,nlayermx) 178 REAL ratiom(ngridmx,nlayermx) 179 180 ! ========================================= 181 165 182 !----------------------------------------------------------------------- 166 183 ! initialisation: … … 172 189 ! zu(:,:)=pu(:,:) 173 190 ! zv(:,:)=pv(:,:) 174 ztv(:,:)=pt(:,:)/zpopsk(:,:) 191 zhc(:,:)=pt(:,:)/zpopsk(:,:) 192 193 ! ********************************************************************** 194 ! Taking into account vertical molar mass gradients 195 ! ********************************************************************** 196 197 if(firstcall) then 198 ico2=0 199 if (tracer) then 200 ! Prepare Special treatment if one of the tracers is CO2 gas 201 do iq=1,nqmx 202 if (noms(iq).eq."co2") then 203 ico2=iq 204 m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol) 205 m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol) 206 ! Compute A and B coefficient use to compute 207 ! mean molecular mass Mair defined by 208 ! 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2 209 ! 1/Mair = A*q(ico2) + B 210 A =(1/m_co2 - 1/m_noco2) 211 B=1/m_noco2 212 end if 213 enddo 214 endif 215 216 firstcall=.false. 217 endif !of if firstcall 218 219 !....................................................................... 220 ! Special treatment for co2 221 !....................................................................... 222 223 if (ico2.ne.0) then 224 ! Special case if one of the tracers is CO2 gas 225 DO l=1,nlayermx 226 DO ig=1,ngridmx 227 ztv(ig,l) = zhc(ig,l)*(A*pq(ig,l,ico2)+B) 228 ENDDO 229 ENDDO 230 else 231 ztv(:,:)=zhc(:,:) 232 end if 233 175 234 176 235 !------------------------------------------------------------------------ … … 332 391 ! couches sont instables. 333 392 !------------------------------------------------------------------------- 334 active (:)=ztv(:,1)>ztv(:,2)393 activecell(:)=ztv(:,1)>ztv(:,2) 335 394 do ig=1,ngridmx 336 395 if (ztv(ig,1)>=(ztv(ig,2))) then … … 380 439 ! Le panache va prendre au debut les caracteristiques de l'air contenu 381 440 ! dans cette couche. 382 if (active (ig)) then441 if (activecell(ig)) then 383 442 ztla(ig,1)=ztv(ig,1) 384 443 !cr: attention, prise en compte de f*(1)=1 => AC : what ? f*(1) =0. ! (d'ou f*(2)=a*(1) … … 407 466 ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test 408 467 do ig=1,ngridmx 409 active (ig)=active(ig) &468 activecell(ig)=activecell(ig) & 410 469 & .and. zw2(ig,l)>1.e-10 & 411 470 & .and. f_star(ig,l)+alim_star(ig,l)>1.e-10 … … 423 482 424 483 do ig=1,ngridmx 425 if(active (ig)) then484 if(activecell(ig)) then 426 485 427 486 ! if(l .lt. lalim(ig)) then … … 453 512 454 513 do ig=1,ngridmx 455 if (active (ig)) then514 if (activecell(ig)) then 456 515 457 516 zw2m=w_est(ig,l+1) … … 459 518 if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then 460 519 entr_star(ig,l)=f_star(ig,l)*zdz* & 461 & MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be) 520 ! & MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be) 521 & MAX(0.,log(1. + 0.027*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1))) 462 522 else 463 523 entr_star(ig,l)=0. … … 480 540 ! new param from continuity eq with a fit on dfdz 481 541 detr_star(ig,l) = f_star(ig,l)*zdz* & 482 & ad 542 ! & ad 543 & Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m) 483 544 484 545 ! & MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005) !svn baseline … … 492 553 else 493 554 detr_star(ig,l)=f_star(ig,l)*zdz* & 494 & bd*zbuoy(ig,l)/zw2m 555 ! & bd*zbuoy(ig,l)/zw2m 556 & Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m) 557 495 558 496 559 ! & MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005) !svn baseline … … 529 592 530 593 DO tic=0,1 ! internal convergence loop 531 activetmp(:)=active (:) .and. f_star(:,l+1)>1.e-10594 activetmp(:)=activecell(:) .and. f_star(:,l+1)>1.e-10 532 595 do ig=1,ngridmx 533 596 if (activetmp(ig)) then … … 895 958 enddo 896 959 897 ! Les "optimisations" du flux sont desactive es : moins de bidouilles960 ! Les "optimisations" du flux sont desactivecelles : moins de bidouilles 898 961 ! je considere que celles ci ne sont pas justifiees ou trop delicates 899 962 ! pour MARS, d'apres des observations LES. … … 1332 1395 1333 1396 !------------------------------------------------------------------ 1397 ! calcul du transport vertical de traceurs 1398 !------------------------------------------------------------------ 1399 1400 ! We only transport co2 tracer because it is coupled to the scheme through theta_m 1401 ! The rest is transported outside the sub-timestep loop 1402 1403 if (ico2.ne.0) then 1404 ! if (nqmx .ne. 0.) then 1405 do l=1,nlayermx 1406 zdzfull(:,l)=zlev(:,l+1)-zlev(:,l) 1407 enddo 1408 1409 modname='tracer' 1410 call thermcell_dqup(ngridmx,nlayermx,ptimestep & 1411 & ,fm,entr,detr, & 1412 & masse,pq(:,:,ico2),pdqadj(:,:,ico2),modname,zdzfull) 1413 ! endif 1414 1415 ! Compute the ratio between theta and theta_m 1416 1417 do l=1,nlayermx 1418 do ig=1,ngridmx 1419 ratiom(ig,l)=1./(A*(pq(ig,l,ico2)+pdqadj(ig,l,ico2)*ptimestep)+B) 1420 enddo 1421 enddo 1422 else 1423 ratiom(:,:)=1. 1424 endif 1425 1426 1427 !------------------------------------------------------------------ 1334 1428 ! incrementation dt 1335 1429 !------------------------------------------------------------------ … … 1337 1431 do l=1,nlayermx 1338 1432 do ig=1,ngridmx 1339 pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l) 1340 ! pdtadj(ig,l)=zdthladj(ig,l)*zpopsk(ig,l) 1341 enddo 1342 enddo 1343 1344 !------------------------------------------------------------------ 1345 ! calcul du transport vertical de traceurs 1346 !------------------------------------------------------------------ 1347 1348 ! if (nqmx .ne. 0.) then 1349 ! modname='tracer' 1350 ! DO iq=1,nqmx 1351 ! call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr, & 1352 ! & masse,pq(:,:,iq),pdqadj(:,:,iq),ztvd,fm_down,ztv,modname,lmax) 1353 ! 1354 ! ENDDO 1355 ! endif 1433 pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)*ratiom(ig,l) 1434 enddo 1435 enddo 1356 1436 1357 1437 !------------------------------------------------------------------ … … 1378 1458 do l=1,nlayermx-1 1379 1459 do ig=1,ngridmx 1380 teta_th_int(ig,l)=0.5*(ztva(ig,l+1)+ztva(ig,l)) 1381 teta_down_int(ig,l) = 0.5*(ztvd(ig,l+1)+ztvd(ig,l)) 1382 teta_env_int(ig,l)=0.5*(ztv(ig,l+1)+ztv(ig,l)) 1460 teta_th_int(ig,l)=0.5*(ztva(ig,l+1)+ztva(ig,l))*ratiom(ig,l) 1461 teta_down_int(ig,l) = 0.5*(ztvd(ig,l+1)+ztvd(ig,l))*ratiom(ig,l) 1462 teta_env_int(ig,l)=0.5*(ztv(ig,l+1)+ztv(ig,l))*ratiom(ig,l) 1383 1463 enddo 1384 1464 enddo … … 1391 1471 do ig=1,ngridmx 1392 1472 heatFlux(ig,l)=fm(ig,l)*(teta_th_int(ig,l)-teta_env_int(ig,l))/(rhobarz(ig,l)) 1393 !buoyancyOut(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)1394 !buoyancyEst(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)1473 buoyancyOut(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) 1474 buoyancyEst(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) 1395 1475 heatFlux_down(ig,l)=fm_down(ig,l)*(teta_down_int(ig,l)-teta_env_int(ig,l))/rhobarz(ig,l) 1396 1476 enddo
Note: See TracChangeset
for help on using the changeset viewer.
