Changeset 2953 for trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
- Timestamp:
- Apr 28, 2023, 2:31:03 PM (19 months ago)
- File:
-
- 1 edited
-
trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F (modified) (47 diffs)
Legend:
- Unmodified
- Added
- Removed
-
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.
