Changeset 2953
- Timestamp:
- Apr 28, 2023, 2:31:03 PM (19 months ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 6 edited
-
README (modified) (1 diff)
-
libf/phymars/co2condens_mod.F (modified) (23 diffs)
-
libf/phymars/physiq_mod.F (modified) (11 diffs)
-
libf/phymars/rocketduststorm_mod.F90 (modified) (4 diffs)
-
libf/phymars/topmons_mod.F90 (modified) (3 diffs)
-
libf/phymars/vdifc_mod.F (modified) (47 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/README
r2951 r2953 4019 4019 too long for fixed fortran format. 4020 4020 Minor fix for picky gfortran compiler in testphys1d (no .eq. for logicals!) 4021 4022 == 28/04/2023 == RV 4023 Correct writing of variables inertiesoil and fluxgeo following r2919 and r2942 4024 4025 == 28/04/2023 == RV 4026 Adapt vdifc, co2condens, rocketduststorm and topmons routines to the subslope parametrisation. -
trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F
r2942 r2953 8 8 CONTAINS 9 9 10 SUBROUTINE co2condens(ngrid,nlayer,nq, ptimestep,10 SUBROUTINE co2condens(ngrid,nlayer,nq,nslope,ptimestep, 11 11 $ pcapcal,pplay,pplev,ptsrf,pt, 12 12 $ pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, … … 34 34 USE vertical_layers_mod, ONLY: ap, bp 35 35 #endif 36 use comslope_mod, ONLY: subslope_dist,def_slope_mean 36 37 IMPLICIT NONE 37 38 c======================================================================= … … 59 60 INTEGER,INTENT(IN) :: nlayer ! number of vertical layers 60 61 INTEGER,INTENT(IN) :: nq ! number of tracers 62 INTEGER,INTENT(IN) :: nslope ! number of subslope 61 63 62 64 REAL,INTENT(IN) :: ptimestep ! physics timestep (s) 63 REAL,INTENT(IN) :: pcapcal(ngrid )65 REAL,INTENT(IN) :: pcapcal(ngrid,nslope) 64 66 REAL,INTENT(IN) :: pplay(ngrid,nlayer) !mid-layer pressure (Pa) 65 67 REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa) 66 REAL,INTENT(IN) :: ptsrf(ngrid ) ! surface temperature (K)68 REAL,INTENT(IN) :: ptsrf(ngrid,nslope) ! surface temperature (K) 67 69 REAL,INTENT(IN) :: pt(ngrid,nlayer) ! atmospheric temperature (K) 68 70 REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential (m2.s-2) … … 73 75 REAL,INTENT(IN) :: pdv(ngrid,nlayer) ! tendency on meridional wind (m/s2) 74 76 ! from previous physical processes 75 REAL,INTENT(IN) :: pdtsrf(ngrid ) ! tendency on surface temperature from77 REAL,INTENT(IN) :: pdtsrf(ngrid,nslope) ! tendency on surface temperature from 76 78 ! previous physical processes (K/s) 77 79 REAL,INTENT(IN) :: pu(ngrid,nlayer) ! zonal wind (m/s) … … 84 86 REAL,INTENT(IN) :: pcondicea_co2microp(ngrid,nlayer)! tendency due to CO2 condensation (kg/kg.s-1) 85 87 86 REAL,INTENT(INOUT) :: piceco2(ngrid ) ! CO2 ice on the surface (kg.m-2)87 REAL,INTENT(INOUT) :: psolaralb(ngrid,2 ) ! albedo of the surface88 REAL,INTENT(INOUT) :: pemisurf(ngrid ) ! emissivity of the surface88 REAL,INTENT(INOUT) :: piceco2(ngrid,nslope) ! CO2 ice on the surface (kg.m-2) 89 REAL,INTENT(INOUT) :: psolaralb(ngrid,2,nslope) ! albedo of the surface 90 REAL,INTENT(INOUT) :: pemisurf(ngrid,nslope) ! emissivity of the surface 89 91 REAL,INTENT(IN) :: rdust(ngrid,nlayer) ! dust effective radius 90 92 91 93 ! tendencies due to CO2 condensation/sublimation: 92 94 REAL,INTENT(OUT) :: pdtc(ngrid,nlayer) ! tendency on temperature (K/s) 93 REAL,INTENT(OUT) :: pdtsrfc(ngrid ) ! tendency on surface temperature (K/s)95 REAL,INTENT(OUT) :: pdtsrfc(ngrid,nslope) ! tendency on surface temperature (K/s) 94 96 REAL,INTENT(OUT) :: pdpsrf(ngrid) ! tendency on surface pressure (Pa/s) 95 97 REAL,INTENT(OUT) :: pduc(ngrid,nlayer) ! tendency on zonal wind (m.s-2) … … 112 114 REAL ztcond (ngrid,nlayer+1) ! CO2 condensation temperature (atm) 113 115 REAL ztcondsol(ngrid) ! CO2 condensation temperature (surface) 114 REAL zdiceco2(ngrid) 116 REAL zdiceco2(ngrid,nslope) 117 REAL zdiceco2_mesh_avg(ngrid) 115 118 REAL zcondicea(ngrid,nlayer) ! condensation rate in layer l (kg/m2/s) 116 REAL zcondices(ngrid) ! condensation rate on the ground (kg/m2/s) 119 REAL zcondices(ngrid,nslope) ! condensation rate on the ground (kg/m2/s) 120 REAL zcondices_mesh_avg(ngrid) ! condensation rate on the ground (kg/m2/s) 117 121 REAL zfallice(ngrid,nlayer+1) ! amount of ice falling from layer l (kg/m2/s) 118 122 REAL condens_layer(ngrid,nlayer) ! co2clouds: condensation rate in layer l (kg/m2/s) … … 122 126 REAL zu(nlayer),zv(nlayer) 123 127 REAL zqc(nlayer,nq),zq1(nlayer) 124 REAL ztsrf(ngrid )128 REAL ztsrf(ngrid,nslope) 125 129 REAL ztc(nlayer), ztm(nlayer+1) 126 130 REAL zum(nlayer+1) , zvm(nlayer+1) … … 129 133 REAL Sm(nlayer),Smq(nlayer,nq),mixmas,qmix 130 134 REAL availco2 131 LOGICAL condsub(ngrid )132 133 real :: emisref(ngrid )135 LOGICAL condsub(ngrid,nslope) 136 137 real :: emisref(ngrid,nslope) 134 138 135 139 REAL zdq_scav(ngrid,nlayer,nq) ! tendency due to scavenging by co2 … … 170 174 REAL masseq(nlayer),wq(nlayer+1) 171 175 INTEGER ifils,iq2 176 177 c Subslope: 178 179 REAL :: emisref_tmp(ngrid) 180 REAL :: alb_tmp(ngrid,2) ! local 181 REAL :: zcondices_tmp(ngrid) ! local 182 REAL :: piceco2_tmp(ngrid) ! local 183 REAL :: pemisurf_tmp(ngrid)! local 184 LOGICAL :: condsub_tmp(ngrid) !local 185 REAL :: zfallice_tmp(ngrid,nlayer+1) 186 REAL :: condens_layer_tmp(ngrid,nlayer) ! co2clouds: condensation rate in layer l (kg/m2/s) 187 INTEGER :: islope 172 188 c---------------------------------------------------------------------- 173 189 … … 231 247 pdqc(1:ngrid,1:nlayer,1:nq) = 0 232 248 233 zcondices(1:ngrid) = 0. 234 pdtsrfc(1:ngrid) = 0. 249 zcondices(1:ngrid,1:nslope) = 0. 250 zcondices_mesh_avg(1:ngrid)=0. 251 pdtsrfc(1:ngrid,1:nslope) = 0. 235 252 pdpsrf(1:ngrid) = 0. 236 condsub(1:ngrid) = .false. 237 zdiceco2(1:ngrid) = 0. 253 condsub(1:ngrid,1:nslope) = .false. 254 zdiceco2(1:ngrid,1:nslope) = 0. 255 zdiceco2_mesh_avg(1:ngrid)=0. 238 256 239 257 zfallheat=0 … … 301 319 pdtc(ig,l)=0. 302 320 IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN 303 condsub(ig )=.true.321 condsub(ig,:)=.true. 304 322 IF (zfallice(ig,l+1).gt.0) then 305 323 zfallheat=zfallice(ig,l+1)* … … 347 365 condens_column(ig) = sum(condens_layer(ig,:)) 348 366 zfallice(ig,1) = zdqssed_co2(ig) 349 piceco2(ig) = piceco2(ig) + zdqssed_co2(ig)*ptimestep 367 DO islope = 1,nslope 368 piceco2(ig,islope) = piceco2(ig,islope) + 369 & zdqssed_co2(ig)*ptimestep * 370 & cos(pi*def_slope_mean(islope)/180.) 371 ENDDO 350 372 ENDDO 351 373 call write_output("co2condens_zfallice", … … 371 393 ztcondsol(ig)= 372 394 & 1./(bcond-acond*log(.01*vmr_co2(ig,1)*pplev(ig,1))) 373 ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep 395 DO islope = 1,nslope 396 ztsrf(ig,islope) = ptsrf(ig,islope) + 397 & pdtsrf(ig,islope)*ptimestep 398 ENDDO 374 399 ENDDO 375 400 … … 388 413 icap=1 389 414 ENDIF 415 416 DO islope = 1,nslope 417 c Need first to put piceco2_slope(ig,islope) in kg/m^2 flat 418 419 piceco2(ig,islope) = piceco2(ig,islope) 420 & /cos(pi*def_slope_mean(islope)/180.) 421 390 422 c 391 423 c Loop on where we have condensation/ sublimation 392 IF ((ztsrf(ig ) .LT. ztcondsol(ig)) .OR. ! ground cond424 IF ((ztsrf(ig,islope) .LT. ztcondsol(ig)) .OR. ! ground cond 393 425 $ (zfallice(ig,1).NE.0.) .OR. ! falling snow 394 $ ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. ! ground sublim. 395 $ ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN 396 condsub(ig) = .true. 426 $ ((ztsrf(ig,islope) .GT. ztcondsol(ig)) .AND. ! ground sublim. 427 $ ((piceco2(ig,islope)+zfallice(ig,1)*ptimestep) 428 $ .NE. 0.))) THEN 429 condsub(ig,islope) = .true. 397 430 398 431 IF (zfallice(ig,1).gt.0 .and. .not. co2clouds) then … … 406 439 c condensation or partial sublimation of CO2 ice 407 440 c """"""""""""""""""""""""""""""""""""""""""""""" 408 zcondices(ig) = pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) 409 & /(latcond*ptimestep) - zfallheat 410 pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep 441 if(ztsrf(ig,islope).LT. ztcondsol(ig)) then 442 c condensation 443 zcondices(ig,islope)=pcapcal(ig,islope) 444 & *(ztcondsol(ig)-ztsrf(ig,islope)) 445 c & /((latcond+cpp*(zt(ig,1)-ztcondsol(ig)))*ptimestep) 446 & /(latcond*ptimestep) 447 & - zfallheat 448 else 449 c sublimation 450 zcondices(ig,islope)=pcapcal(ig,islope) 451 & *(ztcondsol(ig)-ztsrf(ig,islope)) 452 & /(latcond*ptimestep) 453 & - zfallheat 454 endif 455 pdtsrfc(ig,islope) = (ztcondsol(ig) - ztsrf(ig,islope)) 456 & /ptimestep 411 457 #ifdef MESOSCALE 412 458 print*, "not enough CO2 tracer in 1st layer to condense" … … 421 467 availco2 = pq(ig,1,ico2)*((ap(1)-ap(2))+ 422 468 & (bp(1)-bp(2))*(pplev(ig,1)/g - 423 & (zcondices(ig) + zfallice(ig,1))*ptimestep)) 424 if ((zcondices(ig) + condens_layer(ig,1))*ptimestep 469 & (zcondices(ig,islope) + zfallice(ig,1)) 470 & *ptimestep)) 471 if ((zcondices(ig,islope) + condens_layer(ig,1)) 472 & *ptimestep 425 473 & .gt.availco2) then 426 zcondices(ig ) = availco2/ptimestep -474 zcondices(ig,islope) = availco2/ptimestep - 427 475 & condens_layer(ig,1) 428 pdtsrfc(ig ) = (latcond/pcapcal(ig))*429 & (zcondices(ig )+zfallheat)476 pdtsrfc(ig,islope) = (latcond/pcapcal(ig,islope))* 477 & (zcondices(ig,islope)+zfallheat) 430 478 end if 431 479 #else … … 437 485 #endif 438 486 439 c If the entire CO2 ice layer sublimes 487 c If the entire CO2 ice layer sublimes on the slope 440 488 c """""""""""""""""""""""""""""""""""""""""""""""""""" 441 489 c (including what has just condensed in the atmosphere) 442 490 if (co2clouds) then 443 IF((piceco2(ig )/ptimestep).LE.444 & -zcondices(ig ))THEN445 zcondices(ig ) = -piceco2(ig)/ptimestep446 pdtsrfc(ig )=(latcond/pcapcal(ig)) *447 & (zcondices(ig )+zfallheat)491 IF((piceco2(ig,islope)/ptimestep).LE. 492 & -zcondices(ig,islope))THEN 493 zcondices(ig,islope) = -piceco2(ig,islope)/ptimestep 494 pdtsrfc(ig,islope)=(latcond/pcapcal(ig,islope)) * 495 & (zcondices(ig,islope)+zfallheat) 448 496 END IF 449 497 else 450 IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE. 451 & -zcondices(ig))THEN 452 zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1) 453 pdtsrfc(ig)=(latcond/pcapcal(ig)) * 454 & (zcondices(ig)+zfallheat) 498 IF((piceco2(ig,islope)/ptimestep+zfallice(ig,1)).LE. 499 & -zcondices(ig,islope))THEN 500 zcondices(ig,islope) = -piceco2(ig,islope) 501 & /ptimestep - zfallice(ig,1) 502 pdtsrfc(ig,islope)=(latcond/pcapcal(ig,islope)) * 503 & (zcondices(ig,islope)+zfallheat) 455 504 END IF 456 505 end if 457 506 458 c Changing CO2 ice amount and pressure :507 c Changing CO2 ice amount and pressure per slope: 459 508 c """""""""""""""""""""""""""""""""""" 460 zdiceco2(ig ) = zcondices(ig) +zfallice(ig,1)509 zdiceco2(ig,islope) = zcondices(ig,islope)+zfallice(ig,1) 461 510 & + condens_column(ig) 462 511 if (co2clouds) then 463 512 ! add here only direct condensation/sublimation 464 piceco2(ig) = piceco2(ig) + zcondices(ig)*ptimestep 513 piceco2(ig,islope) = piceco2(ig,islope) + 514 & zcondices(ig,islope)*ptimestep 465 515 else 466 516 ! add here also CO2 ice in the atmosphere 467 piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep 517 piceco2(ig,islope) = piceco2(ig,islope) + 518 & zdiceco2(ig,islope)*ptimestep 468 519 end if 469 pdpsrf(ig) = -zdiceco2(ig)*g 520 521 zcondices_mesh_avg(ig) = zcondices_mesh_avg(ig) + 522 & zcondices(ig,islope)* subslope_dist(ig,islope) 523 524 zdiceco2_mesh_avg(ig) = zdiceco2_mesh_avg(ig) + 525 & zdiceco2(ig,islope)* subslope_dist(ig,islope) 526 527 END IF ! if there is condensation/sublimation 528 529 piceco2(ig,islope) = piceco2(ig,islope) 530 & * cos(pi*def_slope_mean(islope)/180.) 531 532 ENDDO !islope 533 534 pdpsrf(ig) = -zdiceco2_mesh_avg(ig)*g 470 535 471 536 IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN … … 480 545 & 'condensing more than total mass', 1) 481 546 ENDIF 482 END IF ! if there is condensation/sublimation 547 483 548 ENDDO ! of DO ig=1,ngrid 549 484 550 485 551 c ******************************************************************** … … 489 555 ! Check that amont of CO2 ice is not problematic 490 556 DO ig=1,ngrid 491 if(.not.piceco2(ig).ge.0.) THEN 492 if(piceco2(ig).le.-5.e-8) print*, 493 $ 'WARNING co2condens piceco2(',ig,')=', piceco2(ig) 494 piceco2(ig)=0. 557 DO islope = 1,nslope 558 if(.not.piceco2(ig,islope).ge.0.) THEN 559 if(piceco2(ig,islope).le.-5.e-8) print*, 560 $ 'WARNING co2condens piceco2(',ig,')=', piceco2(ig,islope) 561 piceco2(ig,islope)=0. 495 562 endif 563 ENDDO 496 564 ENDDO 497 565 498 566 ! Set albedo and emissivity of the surface 499 567 ! ---------------------------------------- 500 CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref) 568 DO islope = 1,nslope 569 piceco2_tmp(:) = piceco2(:,islope) 570 alb_tmp(:,:) = psolaralb(:,:,islope) 571 emisref_tmp(:) = 0. 572 CALL albedocaps(zls,ngrid,piceco2_tmp,alb_tmp,emisref_tmp) 573 psolaralb(:,1,islope) = alb_tmp(:,1) 574 psolaralb(:,2,islope) = alb_tmp(:,2) 575 emisref(:,islope) = emisref_tmp(:) 576 ENDDO 501 577 502 578 ! set pemisurf() to emissiv when there is bare surface (needed for co2snow) 503 579 DO ig=1,ngrid 504 if (piceco2(ig).eq.0) then 505 pemisurf(ig)=emissiv 506 endif 580 DO islope = 1,nslope 581 if (piceco2(ig,islope).eq.0) then 582 pemisurf(ig,islope)=emissiv 583 endif 584 ENDDO 507 585 ENDDO 508 586 … … 515 593 516 594 DO ig=1,ngrid 517 if ( condsub(ig)) then595 if (any(condsub(ig,:))) then 518 596 do l=1,nlayer 519 597 ztc(l) =zt(ig,l) +pdtc(ig,l) *ptimestep … … 527 605 c Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up) 528 606 c """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" 529 zmflux(1) = -zcondices (ig) - zdqssed_co2(ig)607 zmflux(1) = -zcondices_mesh_avg(ig) - zdqssed_co2(ig) 530 608 DO l=1,nlayer 531 609 zmflux(l+1) = zmflux(l) - condens_layer(ig,l) … … 702 780 c CO2 snow / clouds scheme 703 781 c *************************************************************** 704 call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev, 705 & condens_layer,zcondices,zfallice,pemisurf) 782 DO islope = 1,nslope 783 emisref_tmp(:) = emisref(:,islope) 784 condsub_tmp(:) = condsub(:,islope) 785 condens_layer_tmp(:,:) = condens_layer(:,:)* 786 & cos(pi*def_slope_mean(islope)/180.) 787 zcondices_tmp(:) = zcondices(:,islope)* 788 & cos(pi*def_slope_mean(islope)/180.) 789 zfallice_tmp(:,:) = zfallice(:,:)* 790 & cos(pi*def_slope_mean(islope)/180.) 791 pemisurf_tmp(:) = pemisurf(:,islope) 792 793 call co2snow(ngrid,nlayer,ptimestep,emisref_tmp,condsub_tmp, 794 & pplev,condens_layer_tmp,zcondices_tmp,zfallice_tmp, 795 & pemisurf_tmp) 796 pemisurf(:,islope) = pemisurf_tmp(:) 797 798 ENDDO 706 799 c *************************************************************** 707 800 c Ecriture des diagnostiques … … 739 832 & 1./(bcond-acond*log(.01*vmr_co2(ngrid,1)* 740 833 & (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep))) 741 pdtsrfc(ngrid)=(ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep 834 DO islope = 1,nslope 835 pdtsrfc(ngrid,islope)=(ztcondsol(ngrid)- 836 & ztsrf(ngrid,islope))/ptimestep 837 ENDDO ! islope = 1,nslope 742 838 endif 743 839 endif -
trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
r2942 r2953 675 675 c initialize tracers 676 676 c ~~~~~~~~~~~~~~~~~~ 677 678 677 CALL initracer(ngrid,nq,qsurf) 679 678 … … 681 680 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 682 681 CALL surfini(ngrid,qsurf) 683 684 682 CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit) 685 686 683 c initialize soil 687 684 c ~~~~~~~~~~~~~~~ … … 1170 1167 ELSE ! not calling subslope, nslope might be > 1 1171 1168 DO islope = 1,nslope 1172 IF(def_slope_mean(islope).ge.0.) THEN1173 sl_psi = 0. !Northward slope1174 ELSE1175 sl_psi = 180. !Southward slope1176 ENDIF1177 1169 sl_the=abs(def_slope_mean(islope)) 1178 1170 IF (sl_the .gt. 1e-6) THEN 1171 IF(def_slope_mean(islope).ge.0.) THEN 1172 psi_sl(:) = 0. !Northward slope 1173 ELSE 1174 psi_sl(:) = 180. !Southward slope 1175 ENDIF 1179 1176 DO ig=1,ngrid 1180 1177 ztim1=fluxsurf_dn_sw(ig,1,islope) … … 1232 1229 ENDDO 1233 1230 1234 1235 1231 c Net atmospheric radiative heating rate (K.s-1) 1236 1232 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ … … 1254 1250 c Net radiative surface flux (W.m-2) 1255 1251 c ~~~~~~~~~~~~~~~~~~~~~~~~~~ 1252 1256 1253 c 1257 1254 DO ig=1,ngrid … … 1312 1309 c for radiative transfer 1313 1310 & clearatm,icount,zday,zls, 1314 & tsurf ,qsurf_meshavg(:,igcm_co2),igout,1315 & totstormfract,tauscaling,1311 & tsurf_meshavg,qsurf_meshavg(:,igcm_co2), 1312 & igout,totstormfract,tauscaling, 1316 1313 & dust_rad_adjust,IRtoVIScoef, 1314 & albedo_meshavg,emis_meshavg, 1317 1315 c input sub-grid scale cloud 1318 1316 & clearsky,totcloudfrac, … … 1378 1376 & qsurf_meshavg(:,igcm_co2), 1379 1377 & igout,aerosol,tauscaling,dust_rad_adjust, 1380 & IRtoVIScoef,totstormfract,clearatm, 1378 & IRtoVIScoef,albedo_meshavg,emis_meshavg, 1379 & totstormfract,clearatm, 1381 1380 & clearsky,totcloudfrac, 1382 1381 & nohmons, … … 1413 1412 endif ! end of if (dustinjection.gt.0) 1414 1413 1415 1416 1414 c----------------------------------------------------------------------- 1417 1415 c 4. Gravity wave and subgrid scale topography drag : … … 1471 1469 ENDDO 1472 1470 ENDIF 1471 1473 1472 c ---------------------- 1474 1473 … … 2000 1999 c Sedimentation for co2 clouds tracers are inside co2cloud microtimestep 2001 2000 c Zdqssed isn't 2001 2002 2002 call callsedim(ngrid,nlayer,ptimestep, 2003 2003 & zplev,zzlev,zzlay,pt,pdt, … … 2180 2180 zdqc(:,:,:) = 0. 2181 2181 zdqsc(:,:,:) = 0. 2182 CALL co2condens(ngrid,nlayer,nq, ptimestep,2182 CALL co2condens(ngrid,nlayer,nq,nslope,ptimestep, 2183 2183 $ capcal,zplay,zplev,tsurf,pt, 2184 2184 $ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, -
trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90
r2934 r2953 27 27 tsurf,co2ice,igout,totstormfract, & 28 28 tauscaling,dust_rad_adjust, & 29 IRtoVIScoef, 29 IRtoVIScoef,albedo,emis, & 30 30 ! input sub-grid scale cloud 31 31 clearsky,totcloudfrac, & … … 40 40 ,rho_dust 41 41 USE comcstfi_h, only: r,g,cpp,rcp 42 USE dimradmars_mod, only: albedo,naerkind42 USE dimradmars_mod, only: naerkind 43 43 USE comsaison_h, only: dist_sol,mu0,fract 44 USE surfdat_h, only: emis,zmea, zstd, zsig, hmons44 USE surfdat_h, only: zmea, zstd, zsig, hmons 45 45 USE callradite_mod, only: callradite 46 46 use write_output_mod, only: write_output … … 78 78 REAL, INTENT(IN) :: zls 79 79 REAL, INTENT(IN) :: tsurf(ngrid) 80 REAL, INTENT(IN) :: albedo(ngrid,2) 81 REAL, INTENT(IN) :: emis(ngrid) 80 82 REAL,INTENT(IN) :: co2ice(ngrid) ! co2 ice surface layer (kg.m-2) 81 83 INTEGER, INTENT(IN) :: igout … … 255 257 ! 1. Call the second radiative transfer for stormdust, obtain the extra heating 256 258 ! ********************************************************************* 259 257 260 CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo, & 258 261 emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, & -
trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90
r2934 r2953 22 22 icount,zday,zls,tsurf,co2ice,igout, & 23 23 aerosol,tauscaling,dust_rad_adjust, & 24 IRtoVIScoef, 24 IRtoVIScoef,albedo,emis, & 25 25 ! input sub-grid scale rocket dust storm 26 26 totstormfract,clearatm, & … … 37 37 ,rho_dust 38 38 USE comcstfi_h, only: r,g,cpp,rcp 39 USE dimradmars_mod, only: albedo,naerkind39 USE dimradmars_mod, only: naerkind 40 40 USE comsaison_h, only: dist_sol,mu0,fract 41 USE surfdat_h, only: emis,hmons,summit,alpha_hmons, &41 USE surfdat_h, only: hmons,summit,alpha_hmons, & 42 42 hsummit,contains_mons 43 43 USE callradite_mod, only: callradite … … 74 74 REAL, INTENT(IN) :: zls 75 75 REAL, INTENT(IN) :: tsurf(ngrid) 76 REAL, INTENT(IN) :: albedo(ngrid,2) 77 REAL, INTENT(IN) :: emis(ngrid) 76 78 REAL,INTENT(IN) :: co2ice(ngrid) ! co2 ice surface layer (kg.m-2) 77 79 INTEGER, INTENT(IN) :: igout -
trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
r2934 r2953 21 21 & igcm_stormdust_mass, igcm_stormdust_number 22 22 use surfdat_h, only: watercaptag, frost_albedo_threshold, dryness 23 USE comcstfi_h, ONLY: cpp, r, rcp, g 23 USE comcstfi_h, ONLY: cpp, r, rcp, g, pi 24 24 use watersat_mod, only: watersat 25 25 use turb_mod, only: turb_resolved, ustar, tstar … … 29 29 use dust_param_mod, only: doubleq, submicron, lifting 30 30 use write_output_mod, only: write_output 31 use comslope_mod, ONLY: nslope,def_slope,def_slope_mean, 32 & subslope_dist,major_slope,iflat 31 33 32 34 IMPLICIT NONE … … 65 67 REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay) 66 68 REAL,INTENT(IN) :: ph(ngrid,nlay) 67 REAL,INTENT(IN) :: ptsrf(ngrid ),pemis(ngrid)69 REAL,INTENT(IN) :: ptsrf(ngrid,nslope),pemis(ngrid,nslope) 68 70 REAL,INTENT(IN) :: pdufi(ngrid,nlay),pdvfi(ngrid,nlay) 69 71 REAL,INTENT(IN) :: pdhfi(ngrid,nlay) 70 REAL,INTENT(IN) :: pfluxsrf(ngrid )72 REAL,INTENT(IN) :: pfluxsrf(ngrid,nslope) 71 73 REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay) 72 REAL,INTENT(OUT) :: pdtsrf(ngrid ),pdhdif(ngrid,nlay)73 REAL,INTENT(IN) :: pcapcal(ngrid )74 REAL,INTENT(OUT) :: pdtsrf(ngrid,nslope),pdhdif(ngrid,nlay) 75 REAL,INTENT(IN) :: pcapcal(ngrid,nslope) 74 76 REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1) 75 77 … … 87 89 c Traceurs : 88 90 integer,intent(in) :: nq 89 REAL,INTENT(IN) :: pqsurf(ngrid,nq )91 REAL,INTENT(IN) :: pqsurf(ngrid,nq,nslope) 90 92 REAL :: zqsurf(ngrid) ! temporary water tracer 91 93 real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) 92 94 real,intent(out) :: pdqdif(ngrid,nlay,nq) 93 real,intent(out) :: pdqsdif(ngrid,nq )95 real,intent(out) :: pdqsdif(ngrid,nq,nslope) 94 96 REAL,INTENT(in) :: dustliftday(ngrid) 95 97 REAL,INTENT(in) :: local_time(ngrid) … … 100 102 REAL :: pt(ngrid,nlay) 101 103 102 INTEGER ilev,ig,ilay,nlev 104 INTEGER ilev,ig,ilay,nlev,islope 103 105 104 106 REAL z4st,zdplanck(ngrid) … … 106 108 REAL zkq(ngrid,nlay+1) 107 109 REAL zcdv(ngrid),zcdh(ngrid) 108 REAL zcdv_true(ngrid),zcdh_true(ngrid) ! drag coeff are used by the LES to recompute u* and hfx110 REAL, INTENT(IN) :: zcdv_true(ngrid),zcdh_true(ngrid) ! drag coeff are used by the LES to recompute u* and hfx 109 111 REAL zu(ngrid,nlay),zv(ngrid,nlay) 110 112 REAL zh(ngrid,nlay) … … 136 138 REAL zdqsdif(ngrid) ! subtimestep pdqsdif for water ice 137 139 REAL ztsrf(ngrid) ! temporary surface temperature in tsub 138 REAL zdtsrf(ngrid ) ! surface temperature tendancy in tsub139 REAL surf_h2o_lh(ngrid ) ! Surface h2o latent heat flux140 REAL zsurf_h2o_lh(ngrid ) ! Tsub surface h2o latent heat flux140 REAL zdtsrf(ngrid,nslope) ! surface temperature tendancy in tsub 141 REAL surf_h2o_lh(ngrid,nslope) ! Surface h2o latent heat flux 142 REAL zsurf_h2o_lh(ngrid,nslope) ! Tsub surface h2o latent heat flux 141 143 142 144 c For latent heat release from ground water ice sublimation … … 154 156 REAL qsat(ngrid) 155 157 156 REAL hdoflux(ngrid) ! value of vapour flux of HDO 158 REAL hdoflux(ngrid,nslope) ! value of vapour flux of HDO 159 REAL hdoflux_meshavg(ngrid) ! value of vapour flux of HDO 157 160 ! REAL h2oflux(ngrid) ! value of vapour flux of H2O 158 161 REAL old_h2o_vap(ngrid) ! traceur d'eau avant traitement 162 REAL saved_h2o_vap(ngrid) ! traceur d'eau avant traitement 159 163 160 164 REAL kmixmin 161 165 162 166 c Argument added for surface water ice budget: 163 REAL,INTENT(IN) :: watercap(ngrid )164 REAL,INTENT(OUT) :: dwatercap_dif(ngrid )167 REAL,INTENT(IN) :: watercap(ngrid,nslope) 168 REAL,INTENT(OUT) :: dwatercap_dif(ngrid,nslope) 165 169 166 170 c Subtimestep to compute h2o latent heat flux: … … 198 202 !!MARGAUX 199 203 REAL DoH_vap(ngrid,nlay) 204 !! Sub-grid scale slopes 205 REAL :: pdqsdif_tmp(ngrid,nq) ! Temporary for dust lifting 206 REAL :: watercap_tmp(ngrid) 207 REAL :: zq_slope_vap(ngrid,nlay,nq,nslope) 208 REAL :: zq_tmp_vap(ngrid,nlay,nq) 209 REAL :: ptsrf_tmp(ngrid) 210 REAL :: pqsurf_tmp(ngrid,nq) 211 REAL :: pdqsdif_tmphdo(ngrid,nq) 212 REAL :: qsat_tmp(ngrid) 213 INTEGER :: indmax 214 215 character*2 str2 200 216 201 217 c ** un petit test de coherence … … 232 248 ENDIF 233 249 250 DO ig = 1,ngrid 251 ptsrf_tmp(ig) = ptsrf(ig,major_slope(ig)) 252 pqsurf_tmp(ig,:) = pqsurf(ig,:,major_slope(ig)) 253 ENDDO 234 254 235 255 c----------------------------------------------------------------------- … … 243 263 pdvdif(1:ngrid,1:nlay)=0 244 264 pdhdif(1:ngrid,1:nlay)=0 245 pdtsrf(1:ngrid )=0246 zdtsrf(1:ngrid )=0247 surf_h2o_lh(1:ngrid )=0248 zsurf_h2o_lh(1:ngrid )=0265 pdtsrf(1:ngrid,1:nslope)=0 266 zdtsrf(1:ngrid,1:nslope)=0 267 surf_h2o_lh(1:ngrid,1:nslope)=0 268 zsurf_h2o_lh(1:ngrid,1:nslope)=0 249 269 pdqdif(1:ngrid,1:nlay,1:nq)=0 250 pdqsdif(1:ngrid,1:nq)=0 270 pdqsdif(1:ngrid,1:nq,1:nslope)=0 271 pdqsdif_tmp(1:ngrid,1:nq)=0 251 272 zdqsdif(1:ngrid)=0 252 dwatercap_dif(1:ngrid )=0273 dwatercap_dif(1:ngrid,1:nslope)=0 253 274 254 275 c ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp … … 275 296 ENDDO 276 297 DO ig=1,ngrid 277 zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf (ig))298 zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf_tmp(ig)) 278 299 ENDDO 279 300 … … 354 375 ENDDO 355 376 ENDDO 377 356 378 zq(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq)+ 357 379 & pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep … … 366 388 c --------------------- 367 389 368 CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wstar,ptsrf ,ph369 & 390 CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wstar,ptsrf_tmp 391 & ,ph,zcdv_true,zcdh_true) 370 392 371 393 zu2(:)=pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1) … … 383 405 DO ig=1,ngrid 384 406 IF (zcdh_true(ig) .ne. 0.) THEN ! When Cd=Ch=0, u*=t*=0 385 tstar(ig)=(ph(ig,1)-ptsrf(ig))*zcdh(ig)/ustar(ig) 407 tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig)) 408 & *zcdh(ig)/ustar(ig) 386 409 ENDIF 387 410 ENDDO … … 391 414 zcdh(:)=zcdh_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance 392 415 ustar(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)) 393 tstar(:)=(ph(:,1)-ptsrf(:))*zcdh_true(:)/sqrt(zcdv_true(:)) 416 tstar(:)=(ph(:,1)-ptsrf_tmp(:)) 417 & *zcdh_true(:)/sqrt(zcdv_true(:)) 394 418 ENDIF 395 419 … … 454 478 ENDIF 455 479 456 457 458 459 480 c----------------------------------------------------------------------- 460 481 c 4. inversion pour l'implicite sur u … … 498 519 ENDDO 499 520 500 501 502 503 504 521 c----------------------------------------------------------------------- 505 522 c 5. inversion pour l'implicite sur v … … 539 556 ENDDO 540 557 ENDDO 541 542 543 544 545 558 546 559 c----------------------------------------------------------------------- … … 575 588 IF (tke_heat_flux .eq. 0.) THEN 576 589 DO ig=1,ngrid 577 zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig) 590 indmax = major_slope(ig) 591 zdplanck(ig)=z4st*pemis(ig,indmax)*ptsrf(ig,indmax)* 592 & ptsrf(ig,indmax)*ptsrf(ig,indmax) 578 593 ENDDO 579 594 ELSE … … 603 618 ENDIF 604 619 620 621 605 622 DO ig=1,ngrid 606 623 indmax = major_slope(ig) 607 624 ! Initialization of z1 and zd, which do not depend on dmice 608 625 … … 649 666 c ---------------------------------------------------- 650 667 651 z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1) 652 s +zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep 653 z2(ig)= pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig) 668 z1(ig)=pcapcal(ig,indmax)*ptsrf(ig,indmax) 669 s + cpp*zb(ig,1)*zc(ig,1) 670 s + zdplanck(ig)*ptsrf(ig,indmax) 671 s + pfluxsrf(ig,indmax)*ptimestep 672 z2(ig)= pcapcal(ig,indmax)+cpp*zb(ig,1)*(1.-zd(ig,1)) 673 s +zdplanck(ig) 654 674 ztsrf2(ig)=z1(ig)/z2(ig) 655 675 ! pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep !incremented outside loop … … 687 707 ENDDO 688 708 ENDIF 689 709 690 710 ENDDO!of Do j=1,XXX 691 711 pdtsrf(ig,indmax)=(ztsrf2(ig)-ptsrf(ig,indmax))/ptimestep 692 712 ENDDO !of Do ig=1,ngrid 693 694 pdtsrf(:)=(ztsrf2(:)-ptsrf(:))/ptimestep695 713 696 714 DO ig=1,ngrid ! computing sensible heat flux (atm => surface) 697 715 sensibFlux(ig)=cpp*zb(ig,1)/ptimestep*(zhs(ig,1)-ztsrf2(ig)) 698 716 ENDDO 717 718 c Now implicit sheme on each sub-grid subslope: 719 IF (nslope.ne.1) then 720 DO islope=1,nslope 721 DO ig=1,ngrid 722 IF(islope.ne.major_slope(ig)) then 723 IF (tke_heat_flux .eq. 0.) THEN 724 zdplanck(ig)=z4st*pemis(ig,islope)*ptsrf(ig,islope)**3 725 ELSE 726 zdplanck(ig) = 0. 727 ENDIF 728 z1(ig)=pcapcal(ig,islope)*ptsrf(ig,islope) 729 s + cpp*zb(ig,1)*zc(ig,1) 730 s + zdplanck(ig)*ptsrf(ig,islope) 731 s + pfluxsrf(ig,islope)*ptimestep 732 z2(ig)= pcapcal(ig,islope)+cpp*zb(ig,1)*(1.-zd(ig,1)) 733 s +zdplanck(ig) 734 ztsrf2(ig)=z1(ig)/z2(ig) 735 pdtsrf(ig,islope)=(ztsrf2(ig)-ptsrf(ig,islope))/ptimestep 736 ENDIF ! islope != indmax 737 ENDDO ! ig 738 ENDDO !islope 739 ENDIF !nslope.ne.1 699 740 700 741 c----------------------------------------------------------------------- … … 727 768 do ig=1,ngrid 728 769 c if(qsurf(ig,igcm_co2).lt.1) then 729 pdqsdif (ig,igcm_dust_mass) =770 pdqsdif_tmp(ig,igcm_dust_mass) = 730 771 & -alpha_lift(igcm_dust_mass) 731 pdqsdif (ig,igcm_dust_number) =772 pdqsdif_tmp(ig,igcm_dust_number) = 732 773 & -alpha_lift(igcm_dust_number) 733 pdqsdif (ig,igcm_dust_submicron) =774 pdqsdif_tmp(ig,igcm_dust_submicron) = 734 775 & -alpha_lift(igcm_dust_submicron) 735 776 c end if … … 739 780 !or 2 (injection in CL) 740 781 do ig=1,ngrid 741 if(pqsurf (ig,igcm_co2).lt.1) then ! pas de soulevement si glace CO2742 pdqsdif (ig,igcm_dust_mass) =782 if(pqsurf_tmp(ig,igcm_co2).lt.1) then ! pas de soulevement si glace CO2 783 pdqsdif_tmp(ig,igcm_dust_mass) = 743 784 & -alpha_lift(igcm_dust_mass) 744 pdqsdif (ig,igcm_dust_number) =785 pdqsdif_tmp(ig,igcm_dust_number) = 745 786 & -alpha_lift(igcm_dust_number) 746 787 end if … … 748 789 elseif(dustinjection.eq.1)then ! dust injection scheme = 1 injection from surface 749 790 do ig=1,ngrid 750 if(pqsurf (ig,igcm_co2).lt.1) then ! pas de soulevement si glace CO2791 if(pqsurf_tmp(ig,igcm_co2).lt.1) then ! pas de soulevement si glace CO2 751 792 IF((ti_injection_sol.LE.local_time(ig)).and. 752 793 & (local_time(ig).LE.tf_injection_sol)) THEN 753 794 if (rdstorm) then !Rocket dust storm scheme 754 pdqsdif (ig,igcm_stormdust_mass) =795 pdqsdif_tmp(ig,igcm_stormdust_mass) = 755 796 & -alpha_lift(igcm_stormdust_mass) 756 797 & *dustliftday(ig) 757 pdqsdif (ig,igcm_stormdust_number) =798 pdqsdif_tmp(ig,igcm_stormdust_number) = 758 799 & -alpha_lift(igcm_stormdust_number) 759 800 & *dustliftday(ig) 760 pdqsdif (ig,igcm_dust_mass)= 0.761 pdqsdif (ig,igcm_dust_number)= 0.801 pdqsdif_tmp(ig,igcm_dust_mass)= 0. 802 pdqsdif_tmp(ig,igcm_dust_number)= 0. 762 803 else 763 pdqsdif (ig,igcm_dust_mass)=804 pdqsdif_tmp(ig,igcm_dust_mass)= 764 805 & -dustliftday(ig)* 765 806 & alpha_lift(igcm_dust_mass) 766 pdqsdif (ig,igcm_dust_number)=807 pdqsdif_tmp(ig,igcm_dust_number)= 767 808 & -dustliftday(ig)* 768 809 & alpha_lift(igcm_dust_number) 769 810 endif 770 811 if (submicron) then 771 pdqsdif (ig,igcm_dust_submicron) = 0.812 pdqsdif_tmp(ig,igcm_dust_submicron) = 0. 772 813 endif ! if (submicron) 773 814 ELSE ! outside dust injection time frame 774 pdqsdif (ig,igcm_dust_mass)= 0.775 pdqsdif (ig,igcm_dust_number)= 0.815 pdqsdif_tmp(ig,igcm_dust_mass)= 0. 816 pdqsdif_tmp(ig,igcm_dust_number)= 0. 776 817 if (rdstorm) then 777 pdqsdif (ig,igcm_stormdust_mass)= 0.778 pdqsdif (ig,igcm_stormdust_number)= 0.818 pdqsdif_tmp(ig,igcm_stormdust_mass)= 0. 819 pdqsdif_tmp(ig,igcm_stormdust_number)= 0. 779 820 end if 780 821 ENDIF … … 785 826 else if (submicron) then 786 827 do ig=1,ngrid 787 pdqsdif (ig,igcm_dust_submicron) =828 pdqsdif_tmp(ig,igcm_dust_submicron) = 788 829 & -alpha_lift(igcm_dust_submicron) 789 830 end do … … 791 832 #endif 792 833 call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh, 793 & pqsurf (:,igcm_co2),pdqsdif)834 & pqsurf_tmp(:,igcm_co2),pdqsdif_tmp) 794 835 #ifndef MESOSCALE 795 836 endif !doubleq.AND.submicron 796 837 #endif 797 838 else 798 pdqsdif (1:ngrid,1:nq) = 0.839 pdqsdif_tmp(1:ngrid,1:nq) = 0. 799 840 end if 800 841 … … 847 888 zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ 848 889 $ zb(ig,2)*zc(ig,2) + 849 $ (-pdqsdif (ig,iq)) *ptimestep) *z1(ig) !tracer flux from surface890 $ (-pdqsdif_tmp(ig,iq)) *ptimestep) *z1(ig) !tracer flux from surface 850 891 ENDDO 851 892 endif !((iq.eq.igcm_h2o_ice) … … 857 898 ENDDO 858 899 ENDDO 900 DO islope = 1,nslope 901 DO ig = 1,ngrid 902 pdqsdif(ig,iq,islope) = pdqsdif_tmp(ig,iq) 903 & * cos(pi*def_slope_mean(islope)/180.) 904 ENDDO 905 ENDDO 906 859 907 endif! ((.not. water).or.(.not. iq.eq.igcm_h2o_vap)) then 860 908 enddo ! of do iq=1,nq … … 867 915 c de decrire le flux de chaleur latente 868 916 869 870 917 do iq=1,nq 871 918 if ((water).and.(iq.eq.igcm_h2o_vap)) then 872 919 873 920 DO islope = 1,nslope 874 921 DO ig=1,ngrid 875 zqsurf(ig)=pqsurf(ig,igcm_h2o_ice) 922 zqsurf(ig)=pqsurf(ig,igcm_h2o_ice,islope)/ 923 & cos(pi*def_slope_mean(islope)/180.) 924 watercap_tmp(ig) = watercap(ig,islope)/ 925 & cos(pi*def_slope_mean(islope)/180.) 876 926 ENDDO ! ig=1,ngrid 877 927 … … 879 929 c la subroutine est a la fin du fichier 880 930 881 call make_tsub(ngrid,pdtsrf ,zqsurf,931 call make_tsub(ngrid,pdtsrf(:,islope),zqsurf, 882 932 & ptimestep,dtmax,watercaptag, 883 933 & nsubtimestep) … … 886 936 c initialization of ztsrf, which is surface temperature in 887 937 c the subtimestep. 938 saved_h2o_vap(:)= zq(:,1,igcm_h2o_vap) 939 888 940 DO ig=1,ngrid 889 941 subtimestep = ptimestep/nsubtimestep(ig) 890 ztsrf(ig)=ptsrf(ig ) ! +pdtsrf(ig)*subtimestep891 942 ztsrf(ig)=ptsrf(ig,islope) ! +pdtsrf(ig)*subtimestep 943 zq_tmp_vap(ig,:,:) =zq(ig,:,:) 892 944 c Debut du sous pas de temps 893 945 … … 895 947 896 948 c C'est parti ! 897 898 949 zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay) 899 950 & /float(nsubtimestep(ig)) … … 903 954 904 955 z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) 905 zc(ig,nlay)=za(ig,nlay)*zq (ig,nlay,iq)*z1(ig)956 zc(ig,nlay)=za(ig,nlay)*zq_tmp_vap(ig,nlay,iq)*z1(ig) 906 957 zd(ig,nlay)=zb(ig,nlay)*z1(ig) 907 958 DO ilay=nlay-1,2,-1 908 959 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ 909 960 $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) 910 zc(ig,ilay)=(za(ig,ilay)*zq (ig,ilay,iq)+961 zc(ig,ilay)=(za(ig,ilay)*zq_tmp_vap(ig,ilay,iq)+ 911 962 $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) 912 963 zd(ig,ilay)=zb(ig,ilay)*z1(ig) … … 914 965 z1(ig)=1./(za(ig,1)+zb(ig,1)+ 915 966 $ zb(ig,2)*(1.-zd(ig,2))) 916 zc(ig,1)=(za(ig,1)*zq (ig,1,iq)+967 zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+ 917 968 $ zb(ig,2)*zc(ig,2)) * z1(ig) 918 969 919 970 call watersat(1,ztsrf(ig),pplev(ig,1),qsat(ig)) 920 old_h2o_vap(ig)=zq (ig,1,igcm_h2o_vap)971 old_h2o_vap(ig)=zq_tmp_vap(ig,1,igcm_h2o_vap) 921 972 zd(ig,1)=zb(ig,1)*z1(ig) 922 973 zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig) … … 925 976 & *(zq1temp(ig)-qsat(ig)) 926 977 c write(*,*)'subliming more than available frost: qsurf!' 978 927 979 if(.not.watercaptag(ig)) then 928 980 if ((-zdqsdif(ig)*subtimestep) … … 935 987 c write(*,*)'flux vers le sol=',pdqsdif(ig,nq) 936 988 z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2))) 937 zc(ig,1)=(za(ig,1)*zq (ig,1,igcm_h2o_vap)+989 zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,igcm_h2o_vap)+ 938 990 $ zb(ig,2)*zc(ig,2) + 939 991 $ (-zdqsdif(ig)) *subtimestep) *z1(ig) … … 944 996 c Starting upward calculations for water : 945 997 c Actualisation de h2o_vap dans le premier niveau 946 zq (ig,1,igcm_h2o_vap)=zq1temp(ig)998 zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig) 947 999 948 1000 c Take into account the H2O latent heat impact on the surface temperature … … 950 1002 lh=(2834.3-0.28*(ztsrf(ig)-To)- 951 1003 & 0.004*(ztsrf(ig)-To)*(ztsrf(ig)-To))*1.e+3 952 zdtsrf(ig )= zdqsdif(ig)*lh /pcapcal(ig)1004 zdtsrf(ig,islope)= zdqsdif(ig)*lh /pcapcal(ig,islope) 953 1005 endif ! (latentheat_surfwater) then 954 1006 955 1007 DO ilay=2,nlay 956 zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq) 1008 zq_tmp_vap(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay) 1009 & *zq_tmp_vap(ig,ilay-1,iq) 957 1010 ENDDO 958 1011 959 1012 c Subtimestep water budget : 960 1013 961 ztsrf(ig) = ztsrf(ig)+(pdtsrf(ig ) + zdtsrf(ig))962 & 1014 ztsrf(ig) = ztsrf(ig)+(pdtsrf(ig,islope) 1015 & + zdtsrf(ig,islope))*subtimestep 963 1016 zqsurf(ig)= zqsurf(ig)+( 964 1017 & zdqsdif(ig))*subtimestep 965 1018 966 1019 c Monitoring instantaneous latent heat flux in W.m-2 : 967 zsurf_h2o_lh(ig ) = zsurf_h2o_lh(ig)+968 & (zdtsrf(ig)*pcapcal(ig))969 & 1020 zsurf_h2o_lh(ig,islope) = zsurf_h2o_lh(ig,islope)+ 1021 & (zdtsrf(ig,islope)*pcapcal(ig,islope)) 1022 & *subtimestep 970 1023 971 1024 c We ensure that surface temperature can't rise above the solid domain if there … … 973 1026 if(zqsurf(ig) 974 1027 & +zdqsdif(ig)*subtimestep 975 & .gt.frost_albedo_threshold) ! if there is still ice, T cannot exceed To976 & zdtsrf(ig) = min(zdtsrf(ig),(To-ztsrf(ig))/subtimestep) ! ice melt case977 978 1028 & .gt.frost_albedo_threshold) then ! if there is still ice, T cannot exceed To 1029 zdtsrf(ig,islope) = min(zdtsrf(ig,islope), 1030 & (To-ztsrf(ig))/subtimestep) ! ice melt case 1031 endif 979 1032 980 1033 c Fin du sous pas de temps … … 984 1037 c (btw could also compute the post timestep temp and ice 985 1038 c by simply adding the subtimestep trend instead of this) 986 surf_h2o_lh(ig)= zsurf_h2o_lh(ig)/ptimestep 987 pdtsrf(ig)= (ztsrf(ig) - 988 & ptsrf(ig))/ptimestep 989 pdqsdif(ig,igcm_h2o_ice)= 990 & (zqsurf(ig)- pqsurf(ig,igcm_h2o_ice))/ptimestep 1039 surf_h2o_lh(ig,islope)= zsurf_h2o_lh(ig,islope)/ptimestep 1040 pdtsrf(ig,islope)= (ztsrf(ig) - 1041 & ptsrf(ig,islope))/ptimestep 1042 pdqsdif(ig,igcm_h2o_ice,islope)= 1043 & (zqsurf(ig)- pqsurf(ig,igcm_h2o_ice,islope)/ 1044 & cos(pi*def_slope_mean(islope)/180.)) 1045 & /ptimestep 991 1046 992 1047 c if subliming more than qsurf(ice) and on watercaptag, water 993 1048 c sublimates from watercap reservoir (dwatercap_dif is <0) 994 1049 if(watercaptag(ig)) then 995 if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep) 996 & .gt.(pqsurf(ig,igcm_h2o_ice))) then 997 dwatercap_dif(ig)=pdqsdif(ig,igcm_h2o_ice)+ 998 & pqsurf(ig,igcm_h2o_ice)/ptimestep 999 pdqsdif(ig,igcm_h2o_ice)= 1000 & - pqsurf(ig,igcm_h2o_ice)/ptimestep 1050 if ((-pdqsdif(ig,igcm_h2o_ice,islope)*ptimestep) 1051 & .gt.(pqsurf(ig,igcm_h2o_ice,islope) 1052 & /cos(pi*def_slope_mean(islope)/180.))) then 1053 dwatercap_dif(ig,islope)= 1054 & pdqsdif(ig,igcm_h2o_ice,islope)+ 1055 & (pqsurf(ig,igcm_h2o_ice,islope) / 1056 & cos(pi*def_slope_mean(islope)/180.))/ptimestep 1057 pdqsdif(ig,igcm_h2o_ice,islope)= 1058 & - (pqsurf(ig,igcm_h2o_ice,islope)/ 1059 & cos(pi*def_slope_mean(islope)/180.))/ptimestep 1001 1060 endif! ((-pdqsdif(ig)*ptimestep) 1002 1061 endif !(watercaptag(ig)) then 1003 1062 zq_slope_vap(ig,:,:,islope) = zq_tmp_vap(ig,:,:) 1004 1063 ENDDO ! of DO ig=1,ngrid 1064 ENDDO ! islope 1065 1066 c Some grid box averages: interface with the atmosphere 1067 DO ig = 1,ngrid 1068 DO ilay = 1,nlay 1069 zq(ig,ilay,iq) = 0. 1070 DO islope = 1,nslope 1071 zq(ig,ilay,iq) = zq(ig,ilay,iq) + 1072 $ zq_slope_vap(ig,ilay,iq,islope) * 1073 $ subslope_dist(ig,islope) 1074 ENDDO 1075 ENDDO 1076 ENDDO 1077 1078 ! Recompute values in kg/m^2 slopped 1079 DO ig = 1,ngrid 1080 DO islope = 1,nslope 1081 pdqsdif(ig,igcm_h2o_ice,islope) = 1082 & pdqsdif(ig,igcm_h2o_ice,islope) 1083 & * cos(pi*def_slope_mean(islope)/180.) 1084 1085 dwatercap_dif(ig,islope) = 1086 & dwatercap_dif(ig,islope) 1087 & * cos(pi*def_slope_mean(islope)/180.) 1088 ENDDO 1089 ENDDO 1090 1005 1091 END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap)) 1006 1092 … … 1031 1117 ENDDO 1032 1118 ENDDO 1119 hdoflux_meshavg(:) = 0. 1120 DO islope = 1,nslope 1121 1122 pdqsdif_tmphdo(:,:) = pdqsdif(:,:,islope) 1123 & /cos(pi*def_slope_mean(islope)/180.) 1124 1125 call watersat(ngrid,pdtsrf(:,islope)*ptimestep + 1126 & ptsrf(:,islope),pplev(:,1),qsat_tmp) 1033 1127 1034 1128 CALL hdo_surfex(ngrid,nlay,nq,ptimestep, 1035 & zt,pplay,zq,pqsurf, 1036 & old_h2o_vap,qsat,pdqsdif,dwatercap_dif, 1037 & hdoflux) 1129 & zt,pplay,zq,pqsurf(:,:,islope), 1130 & saved_h2o_vap,qsat_tmp, 1131 & pdqsdif_tmphdo, 1132 & dwatercap_dif(:,islope)/cos(pi*def_slope_mean(islope)/180.), 1133 & hdoflux(:,islope)) 1134 1135 pdqsdif(:,:,islope) = pdqsdif_tmphdo(:,:) * 1136 & cos(pi*def_slope_mean(islope)/180.) 1137 DO ig = 1,ngrid 1138 hdoflux_meshavg(ig) = hdoflux_meshavg(ig) + 1139 & hdoflux(ig,islope)*subslope_dist(ig,islope) 1140 1141 ENDDO !ig = 1,ngrid 1142 ENDDO !islope = 1,nslope 1143 1038 1144 DO ig=1,ngrid 1039 1145 z1(ig)=1./(za(ig,1)+zb(ig,1)+ … … 1041 1147 zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ 1042 1148 $ zb(ig,2)*zc(ig,2) + 1043 $ (-hdoflux (ig)) *ptimestep) *z1(ig) !tracer flux from surface1149 $ (-hdoflux_meshavg(ig)) *ptimestep) *z1(ig) !tracer flux from surface 1044 1150 ENDDO 1045 1151 … … 1060 1166 call write_output("surf_h2o_lh", 1061 1167 & "Ground ice latent heat flux", 1062 & "W.m-2",surf_h2o_lh(: ))1168 & "W.m-2",surf_h2o_lh(:,iflat)) 1063 1169 1064 1170 C Diagnostic output for HDO … … 1066 1172 ! CALL write_output('hdoflux', 1067 1173 ! & 'hdoflux', 1068 ! & ' ',hdoflux (:))1174 ! & ' ',hdoflux_meshavg(:)) 1069 1175 ! CALL write_output('h2oflux', 1070 1176 ! & 'h2oflux', … … 1101 1207 WRITE(*,'(a10,3a15)') 1102 1208 s 'theta(t)','theta(t+dt)','u(t)','u(t+dt)' 1103 PRINT*,ptsrf(ngrid/2+1 ),ztsrf2(ngrid/2+1)1209 PRINT*,ptsrf(ngrid/2+1,:),ztsrf2(ngrid/2+1) 1104 1210 DO ilev=1,nlay 1105 1211 WRITE(*,'(4f15.7)')
Note: See TracChangeset
for help on using the changeset viewer.
