Changeset 2515 for trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
- Timestamp:
- May 6, 2021, 12:01:28 PM (4 years ago)
- File:
-
- 1 edited
-
trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
r2409 r2515 26 26 use compute_dtau_mod, only: ti_injection_sol,tf_injection_sol 27 27 use hdo_surfex_mod, only: hdo_surfex 28 c use geometry_mod, only: longitude_deg,latitude_deg ! Joseph 28 29 use dust_param_mod, only: doubleq, submicron, lifting 29 30 … … 86 87 integer,intent(in) :: nq 87 88 REAL,INTENT(IN) :: pqsurf(ngrid,nq) 89 REAL :: zqsurf(ngrid) ! temporary water tracer 88 90 real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) 89 91 real,intent(out) :: pdqdif(ngrid,nlay,nq) … … 126 128 REAL,SAVE :: acond,bcond 127 129 130 c Subtimestep & implicit treatment of water vapor 131 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 132 REAL zdqsdif(ngrid) ! subtimestep pdqsdif for water ice 133 REAL zwatercap(ngrid) ! subtimestep watercap for water ice 134 REAL ztsrf(ngrid) ! temporary surface temperature in tsub 135 REAL zdtsrf(ngrid) ! surface temperature tendancy in tsub 136 128 137 c For latent heat release from ground water ice sublimation 129 REAL tsrf_lh(ngrid) ! temporary surface temperature 138 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 139 REAL tsrf_lh(ngrid) ! temporary surface temperature with lh effect 130 140 REAL lh ! latent heat, formulation given in the Technical Document: 131 141 ! "Modeling water ice sublimation under Phoenix-like conditions", Montmessin et al. 2004 132 ! REAL tsrf_lw(ngrid)133 ! REAL alpha134 ! REAL T1,T2135 ! SAVE T1,T2136 ! DATA T1,T2/-604.3,1080.7/ ! zeros of latent heat equation for ice137 142 138 143 c Tracers : … … 153 158 REAL,INTENT(IN) :: watercap(ngrid) 154 159 REAL,INTENT(OUT) :: dwatercap_dif(ngrid) 160 REAL :: zdwatercap_dif(ngrid) 161 162 c Subtimestep to compute h2o latent heat flux: 163 REAL :: dtmax = 0.5 ! subtimestep temp criterion 164 INTEGER tsub ! adaptative subtimestep (seconds) 165 REAL subtimestep !ptimestep/nsubtimestep 166 INTEGER nsubtimestep(ngrid) ! number of subtimestep (int) 155 167 156 168 c Mass-variation scheme : … … 226 238 pdhdif(1:ngrid,1:nlay)=0 227 239 pdtsrf(1:ngrid)=0 240 zdtsrf(1:ngrid)=0 228 241 pdqdif(1:ngrid,1:nlay,1:nq)=0 229 242 pdqsdif(1:ngrid,1:nq)=0 243 zdqsdif(1:ngrid)=0 244 zwatercap(1:ngrid)=0 230 245 dwatercap_dif(1:ngrid)=0 246 zdwatercap_dif(1:ngrid)=0 231 247 232 248 c ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp … … 788 804 c OU calcul de la valeur de q a la surface (water) : 789 805 c ---------------------------------------- 790 if (water) then791 call watersat(ngrid,ptsrf,pplev(1,1),qsat)792 end if793 806 794 807 c Inversion pour l'implicite sur q 808 c Cas des traceurs qui ne sont pas h2o_vap 809 c h2o_vap est traite plus loin avec un sous pas de temps 810 c hdo_vap est traite ensuite car dependant de h2o_vap 795 811 c -------------------------------- 796 do iq=1,nq !for all tracers including stormdust 812 813 do iq=1,nq !for all tracers except water vapor 814 if ((.not. water).or.(.not. iq.eq.igcm_h2o_vap).or. 815 & (.not. iq.eq.igcm_hdo_vap)) then 816 817 797 818 zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay) 798 799 if ((water).and.(iq.eq.igcm_h2o_vap)) then 800 c This line is required to account for turbulent transport 801 c from surface (e.g. ice) to mid-layer of atmosphere: 802 zb(1:ngrid,1)=zcdv(1:ngrid)*zb0(1:ngrid,1) 803 zb(1:ngrid,1)=dryness(1:ngrid)*zb(1:ngrid,1) 804 else ! (re)-initialize zb(:,1) 805 zb(1:ngrid,1)=0 806 end if 819 zb(1:ngrid,1)=0 807 820 808 821 DO ig=1,ngrid … … 822 835 ENDDO 823 836 824 if ( (water.and.(iq.eq.igcm_h2o_ice)) 825 $ .or. (hdo.and.(iq.eq.igcm_hdo_ice)) ) then 826 ! special case for water ice tracer: do not include 827 ! h2o ice tracer from surface (which is set when handling 828 ! h2o vapour case (see further down). 829 DO ig=1,ngrid 830 z1(ig)=1./(za(ig,1)+zb(ig,1)+ 831 $ zb(ig,2)*(1.-zd(ig,2))) 832 zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ 833 $ zb(ig,2)*zc(ig,2)) *z1(ig) 834 ENDDO 835 else if (hdo.and.(iq.eq.igcm_hdo_vap)) then 837 if ((iq.eq.igcm_h2o_ice) 838 $ .or. (hdo.and.(iq.eq.igcm_hdo_ice) )) then 839 840 DO ig=1,ngrid 841 z1(ig)=1./(za(ig,1)+zb(ig,1)+ 842 $ zb(ig,2)*(1.-zd(ig,2))) 843 zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ 844 $ zb(ig,2)*zc(ig,2)) *z1(ig) !special case h2o_ice 845 ENDDO 846 else ! every other tracer 847 DO ig=1,ngrid 848 z1(ig)=1./(za(ig,1)+zb(ig,1)+ 849 $ zb(ig,2)*(1.-zd(ig,2))) 850 zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ 851 $ zb(ig,2)*zc(ig,2) + 852 $ (-pdqsdif(ig,iq)) *ptimestep) *z1(ig) !tracer flux from surface 853 ENDDO 854 endif !((iq.eq.igcm_h2o_ice) 855 c Starting upward calculations for simple mixing of tracer (dust) 856 DO ig=1,ngrid 857 zq(ig,1,iq)=zc(ig,1) 858 DO ilay=2,nlay 859 zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq) 860 ENDDO 861 ENDDO 862 endif! ((.not. water).or.(.not. iq.eq.igcm_h2o_vap)) then 863 enddo ! of do iq=1,nq 864 865 c --------- h2o_vap -------------------------------- 866 867 868 c Traitement de la vapeur d'eau h2o_vap 869 c Utilisation d'un sous pas de temps afin 870 c de decrire le flux de chaleur latente 871 872 873 do iq=1,nq 874 if ((water).and.(iq.eq.igcm_h2o_vap)) then 875 876 877 DO ig=1,ngrid 878 zqsurf(ig)=pqsurf(ig,igcm_h2o_ice) 879 ENDDO ! ig=1,ngrid 880 881 c make_tsub : sous pas de temps adaptatif 882 c la subroutine est a la fin du fichier 883 884 call make_tsub(ngrid,pdtsrf,zqsurf, 885 & ptimestep,dtmax,watercaptag, 886 & nsubtimestep) 887 888 c Calculation for turbulent exchange with the surface (for ice) 889 c initialization of ztsrf, which is surface temperature in 890 c the subtimestep. 891 DO ig=1,ngrid 892 subtimestep = ptimestep/nsubtimestep(ig) 893 ztsrf(ig)=ptsrf(ig) ! +pdtsrf(ig)*subtimestep 894 zwatercap(ig)=watercap(ig) 895 896 c Debut du sous pas de temps 897 898 DO tsub=1,nsubtimestep(ig) 899 900 c C'est parti ! 901 902 zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay) 903 & /float(nsubtimestep(ig)) 904 zb(1:ngrid,1)=zcdv(1:ngrid)*zb0(1:ngrid,1) 905 & /float(nsubtimestep(ig)) 906 zb(1:ngrid,1)=dryness(1:ngrid)*zb(1:ngrid,1) 907 908 z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) 909 zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig) 910 zd(ig,nlay)=zb(ig,nlay)*z1(ig) 911 DO ilay=nlay-1,2,-1 912 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ 913 $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) 914 zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+ 915 $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) 916 zd(ig,ilay)=zb(ig,ilay)*z1(ig) 917 ENDDO 918 z1(ig)=1./(za(ig,1)+zb(ig,1)+ 919 $ zb(ig,2)*(1.-zd(ig,2))) 920 zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ 921 $ zb(ig,2)*zc(ig,2)) * z1(ig) 922 923 call watersat(1,ztsrf(ig),pplev(1,1),qsat(ig)) 924 old_h2o_vap(ig)=zq(ig,1,igcm_h2o_vap) 925 zd(ig,1)=zb(ig,1)*z1(ig) 926 zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig) 927 928 zdqsdif(ig)=rho(ig)*dryness(ig)*zcdv(ig) 929 & *(zq1temp(ig)-qsat(ig)) 930 if ((-zdqsdif(ig)*subtimestep) 931 & .gt.(zqsurf(ig))) then 932 c write(*,*)'subliming more than available frost: qsurf!' 933 if(watercaptag(ig)) then 934 c dwatercap_dif same sign as pdqsdif (pos. toward ground) 935 zdwatercap_dif(ig) = zdqsdif(ig) 936 & + zqsurf(ig)/subtimestep 937 zdqsdif(ig)= 938 & -zqsurf(ig)/subtimestep 939 c write(*,*) "subliming more than available frost: qsurf!" 940 else ! (case not.watercaptag(ig)) 941 zdqsdif(ig)= 942 & -zqsurf(ig)/subtimestep 943 zdwatercap_dif(ig) = 0.0 944 945 c write(*,*)'flux vers le sol=',pdqsdif(ig,nq) 946 z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2))) 947 zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+ 948 $ zb(ig,2)*zc(ig,2) + 949 $ (-zdqsdif(ig)-zdwatercap_dif(ig)) *subtimestep) *z1(ig) 950 zq1temp(ig)=zc(ig,1) 951 endif !if .not.watercaptag(ig) 952 endif ! if sublim more than surface 953 954 c Starting upward calculations for water : 955 c Actualisation de h2o_vap dans le premier niveau 956 zq(ig,1,igcm_h2o_vap)=zq1temp(ig) 957 958 c Take into account the H2O latent heat impact on the surface temperature 959 if (latentheat_surfwater) then 960 lh=(2834.3-0.28*(ztsrf(ig)-To)- 961 & 0.004*(ztsrf(ig)-To)*(ztsrf(ig)-To))*1.e+3 962 zdtsrf(ig)= (zdwatercap_dif(ig)+ 963 & zdqsdif(ig))*lh /pcapcal(ig) 964 endif ! (latentheat_surfwater) then 965 966 967 c Subtimestep water budget : 968 969 ztsrf(ig) = ztsrf(ig)+(pdtsrf(ig) + zdtsrf(ig)) 970 & *subtimestep 971 zqsurf(ig)= zqsurf(ig)+( 972 & zdqsdif(ig))*subtimestep 973 zwatercap(ig)= zwatercap(ig)+ 974 & zdwatercap_dif(ig)*subtimestep 975 976 977 c We ensure that surface temperature can't rise above the solid domain if there 978 c is still ice on the surface (oldschool) 979 if(zqsurf(ig) 980 & +zdqsdif(ig)*subtimestep 981 & .gt.frost_albedo_threshold) ! if there is still ice, T cannot exceed To 982 & zdtsrf(ig) = min(zdtsrf(ig),(To-ztsrf(ig))/subtimestep) ! ice melt case 983 984 985 DO ilay=2,nlay 986 zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq) 987 ENDDO 988 989 c Fin du sous pas de temps 990 991 ENDDO ! tsub=1,nsubtimestep 992 993 c Integration of subtimestep water budget : 994 995 pdtsrf(ig)= (ztsrf(ig) - 996 & ptsrf(ig))/ptimestep 997 pdqsdif(ig,igcm_h2o_ice)= 998 & (zqsurf(ig)- pqsurf(ig,igcm_h2o_ice))/ptimestep 999 dwatercap_dif(ig)=(zwatercap(ig) - 1000 & watercap(ig))/ptimestep 1001 1002 ENDDO ! of DO ig=1,ngrid 1003 END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap)) 1004 1005 c --------- end of h2o_vap ---------------------------- 1006 1007 c --------- hdo_vap ----------------------------------- 1008 1009 c hdo_ice has already been with along h2o_ice 1010 c amongst "normal" tracers (ie not h2o_vap) 1011 1012 if (hdo.and.(iq.eq.igcm_hdo_vap)) then 1013 zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay) 1014 zb(1:ngrid,1)=0 1015 1016 DO ig=1,ngrid 1017 z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) 1018 zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig) 1019 zd(ig,nlay)=zb(ig,nlay)*z1(ig) 1020 ENDDO 1021 1022 DO ilay=nlay-1,2,-1 1023 DO ig=1,ngrid 1024 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ 1025 $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) 1026 zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+ 1027 $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) 1028 zd(ig,ilay)=zb(ig,ilay)*z1(ig) 1029 ENDDO 1030 ENDDO 1031 836 1032 CALL hdo_surfex(ngrid,nlay,nq,ptimestep, 837 1033 & zt,zq,pqsurf, … … 846 1042 ENDDO 847 1043 848 else ! general case849 1044 DO ig=1,ngrid 850 z1(ig)=1./(za(ig,1)+zb(ig,1)+ 851 $ zb(ig,2)*(1.-zd(ig,2))) 852 zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ 853 $ zb(ig,2)*zc(ig,2) + 854 $ (-pdqsdif(ig,iq)) *ptimestep) *z1(ig) !tracer flux from surface 1045 zq(ig,1,iq)=zc(ig,1) 1046 DO ilay=2,nlay 1047 zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq) 1048 ENDDO 855 1049 ENDDO 856 857 endif ! of if (water.and.(iq.eq.igcm_h2o_ice)) 858 859 IF ((water).and.(iq.eq.igcm_h2o_vap)) then 860 c Calculation for turbulent exchange with the surface (for ice) 861 DO ig=1,ngrid 862 old_h2o_vap(ig)=zq(ig,1,igcm_h2o_vap) 863 zd(ig,1)=zb(ig,1)*z1(ig) 864 zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig) 865 866 pdqsdif(ig,igcm_h2o_ice)=rho(ig)*dryness(ig)*zcdv(ig) 867 & *(zq1temp(ig)-qsat(ig)) 868 c write(*,*)'flux vers le sol=',pdqsdif(ig,nq) 869 END DO 870 871 DO ig=1,ngrid 872 IF (hdo) then !if hdo, need to treat polar caps differently 873 h2oflux(ig) = pdqsdif(ig,igcm_h2o_ice) ! h2oflux is 874 ! uncorrected waterflux 875 ENDIF !hdo 876 877 if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep) 878 & .gt.(pqsurf(ig,igcm_h2o_ice))) then 879 c write(*,*)'subliming more than available frost: qsurf!' 880 if(watercaptag(ig)) then 881 c dwatercap_dif same sign as pdqsdif (pos. toward ground) 882 dwatercap_dif(ig) = pdqsdif(ig,igcm_h2o_ice) 883 & + pqsurf(ig,igcm_h2o_ice)/ptimestep 884 pdqsdif(ig,igcm_h2o_ice)= 885 & -pqsurf(ig,igcm_h2o_ice)/ptimestep 886 c write(*,*) "subliming more than available frost: qsurf!" 887 else ! (case not.watercaptag(ig)) 888 pdqsdif(ig,igcm_h2o_ice)= 889 & -pqsurf(ig,igcm_h2o_ice)/ptimestep 890 dwatercap_dif(ig) = 0.0 891 if (hdo) then 892 h2oflux(ig) = pdqsdif(ig,igcm_h2o_ice) 893 endif 894 895 c write(*,*)'flux vers le sol=',pdqsdif(ig,nq) 896 z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2))) 897 zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+ 898 $ zb(ig,2)*zc(ig,2) + 899 $ (-pdqsdif(ig,igcm_h2o_ice)) *ptimestep) *z1(ig) 900 zq1temp(ig)=zc(ig,1) 901 endif !if .not.watercaptag(ig) 902 endif ! if sublim more than surface 903 904 c Starting upward calculations for water : 905 zq(ig,1,igcm_h2o_vap)=zq1temp(ig) 906 907 !c Take into account H2O latent heat in surface energy budget 908 !c We solve dT/dt = (2834.3-0.28*(T-To)-0.004*(T-To)^2)*1e3*iceflux/cpp 909 ! tsrf_lw(ig) = ptsrf(ig) + pdtsrf(ig) *ptimestep 910 ! 911 ! alpha = exp(-4*abs(T1-T2)*pdqsdif(ig,igcm_h2o_ice) 912 ! & *ptimestep/pcapcal(ig)) 913 ! 914 ! tsrf_lw(ig) = (tsrf_lw(ig)*(T2-alpha*T1)+T1*T2*(alpha-1)) 915 ! & /(tsrf_lw(ig)*(1-alpha)+alpha*T2-T1) ! surface temperature at t+1 916 ! 917 ! pdtsrf(ig) = (tsrf_lw(ig)-ptsrf(ig))/ptimestep 918 c Take into account the H2O latent heat impact on the surface temperature 919 if (latentheat_surfwater) then 920 tsrf_lh(ig) = ptsrf(ig) + pdtsrf(ig) *ptimestep 921 lh=(2834.3-0.28*(tsrf_lh(ig)-To)- 922 & 0.004*(tsrf_lh(ig)-To)*(tsrf_lh(ig)-To))*1.e+3 923 pdtsrf(ig)= pdtsrf(ig) + (dwatercap_dif(ig)+ 924 & pdqsdif(ig,igcm_h2o_ice))*lh /pcapcal(ig) 925 endif 926 927 if(pqsurf(ig,igcm_h2o_ice) 928 & +pdqsdif(ig,igcm_h2o_ice)*ptimestep 929 & .gt.frost_albedo_threshold) ! if there is still ice, T cannot exceed To 930 & pdtsrf(ig) = min(pdtsrf(ig),(To-ptsrf(ig))/ptimestep) ! ice melt case 931 932 ENDDO ! of DO ig=1,ngrid 933 ELSE 934 c Starting upward calculations for simple mixing of tracer (dust) 935 DO ig=1,ngrid 936 zq(ig,1,iq)=zc(ig,1) 937 ENDDO 938 END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap)) 939 940 DO ilay=2,nlay 941 DO ig=1,ngrid 942 zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq) 943 ENDDO 944 ENDDO 1050 endif ! (hdo.and.(iq.eq.igcm_hdo_vap)) 1051 1052 c --------- end of hdo ---------------------------- 1053 945 1054 enddo ! of do iq=1,nq 946 1055 end if ! of if(tracer) 1056 1057 c --------- end of tracers ---------------------------- 1058 947 1059 948 1060 C Diagnostic output for HDO … … 1004 1116 c==================================== 1005 1117 1006 1118 SUBROUTINE make_tsub(naersize,dtsurf,qsurf,ptimestep, 1119 $ dtmax,watercaptag,ntsub) 1120 1121 c Pas de temps adaptatif en estimant le taux de sublimation 1122 c et en adaptant avec un critere "dtmax" du chauffage a accomoder 1123 c dtmax est regle empiriquement (pour l'instant) a 0.5 K 1124 1125 integer,intent(in) :: naersize 1126 real,intent(in) :: dtsurf(naersize) 1127 real,intent(in) :: qsurf(naersize) 1128 logical,intent(in) :: watercaptag(naersize) 1129 real,intent(in) :: ptimestep 1130 real,intent(in) :: dtmax 1131 real :: ztsub(naersize) 1132 integer :: i 1133 integer,intent(out) :: ntsub(naersize) 1134 1135 do i=1,naersize 1136 if ((qsurf(i).eq.0).and. 1137 & (.not.watercaptag(i))) then 1138 ntsub(i) = 1 1139 else 1140 ztsub(i) = ptimestep * dtsurf(i) / dtmax 1141 ntsub(i) = ceiling(abs(ztsub(i))) 1142 endif ! (qsurf(i).eq.0) then 1143 c 1144 c write(78,*), dtsurf*ptimestep, dtsurf, ntsub 1145 enddo! 1=1,ngrid 1146 1147 1148 1149 END SUBROUTINE make_tsub 1007 1150 END MODULE vdifc_mod
Note: See TracChangeset
for help on using the changeset viewer.
