Changeset 3275 for trunk/LMDZ.PLUTO/libf
- Timestamp:
- Mar 20, 2024, 3:05:14 PM (9 months ago)
- Location:
- trunk/LMDZ.PLUTO/libf
- Files:
-
- 2 added
- 19 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/lect_start_archive.F
r3198 r3275 776 776 DO iq=1,nqtot 777 777 txt=trim(tname(iq))//"_surf" 778 if (txt.eq."h2o_ vap") then779 ! There is no surface tracer for h2o_ vap;778 if (txt.eq."h2o_gas") then 779 ! There is no surface tracer for h2o_gas; 780 780 ! "h2o_ice" should be loaded instead 781 781 txt="h2o_ice_surf" 782 782 write(*,*) 'lect_start_archive: loading surface tracer', 783 & ' h2o_ice instead of h2o_ vap'783 & ' h2o_ice instead of h2o_gas' 784 784 endif 785 785 -
trunk/LMDZ.PLUTO/libf/phypluto
-
Property
svn:ignore
set to
*.mod
vdifc_mod.F.save
-
Property
svn:ignore
set to
-
trunk/LMDZ.PLUTO/libf/phypluto/aeropacity.F90
r3195 r3275 94 94 95 95 real CLFtot 96 integer igen_ice,igen_ vap! to store the index of generic tracer96 integer igen_ice,igen_gas ! to store the index of generic tracer 97 97 logical dummy_bool ! dummy boolean just in case we need one 98 98 ! integer i_rgcs_ice(aerogeneric) -
trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90
r3197 r3275 7 7 subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf, & 8 8 albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, & 9 tsurf,fract,dist_star,aerosol,muvar, &9 zzlay,tsurf,fract,dist_star,aerosol,muvar, & 10 10 dtlw,dtsw,fluxsurf_lw, & 11 11 fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw, & … … 30 30 use aeropacity_mod, only: aeropacity 31 31 use aeroptproperties_mod, only: aeroptproperties 32 use tracer_h, only: constants_epsi_generic 33 use comcstfi_mod, only: pi, mugaz, cpp 32 use tracer_h, only: constants_epsi_generic,igcm_ch4_gas,igcm_n2,mmol 33 use comcstfi_mod, only: pi, mugaz, cpp, r, g 34 34 use callkeys_mod, only: varactive,diurnal,tracer,varfixed,satval, & 35 35 diagdtau,kastprof,strictboundcorrk,specOLR, & 36 36 tplanckmin,tplanckmax,global1d, & 37 generic_condensation, aerohaze, haze_radproffix 37 generic_condensation,aerohaze,haze_radproffix,& 38 methane,carbox,cooling,nlte,strobel,& 39 ch4fix,vmrch4_proffix,vmrch4fix 38 40 use optcv_mod, only: optcv 39 41 use optci_mod, only: optci … … 44 46 use generic_tracer_index_mod, only: generic_tracer_index 45 47 use planetwide_mod, only: planetwide_maxval, planetwide_minval 48 use radcommon_h, only: wavev,wavei 46 49 implicit none 47 50 … … 80 83 REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa). 81 84 REAL,INTENT(IN) :: pt(ngrid,nlayer) ! Air temperature (K). 85 REAL,INTENT(IN) :: zzlay(ngrid,nlayer) ! Mid-layer altitude 82 86 REAL,INTENT(IN) :: tsurf(ngrid) ! Surface temperature (K). 83 87 REAL,INTENT(IN) :: fract(ngrid) ! Fraction of day. … … 157 161 !$OMP THREADPRIVATE(tauaero) 158 162 REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn 163 real*8 nfluxtopv_nu(L_NSPECTV) 159 164 REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2). 160 165 REAL*8 nfluxtopi_nu(L_NSPECTI) ! Net band-resolved IR flux at TOA (W/m2). 161 166 REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! For 1D diagnostic. 162 167 REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD) 168 real*8 fmneti_nu(L_NLAYRAD,L_NSPECTI) ! 169 real*8 fmnetv_nu(L_NLAYRAD,L_NSPECTV) ! 163 170 REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD) 164 171 REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD) … … 211 218 REAL*8 surface_stellar_flux ! Stellar flux reaching the surface. Useful for equivalent albedo calculation. 212 219 220 ! NLTE factor for CH4 221 real eps_nlte_sw23(ngrid,nlayer) ! CH4 NLTE efficiency factor for zdtsw 222 real eps_nlte_sw33(ngrid,nlayer) ! CH4 NLTE efficiency factor for zdtsw 223 real eps_nlte_lw(ngrid,nlayer) ! CH4 NLTE efficiency factor for zdtsw 224 integer Nfine,ifine 225 parameter(Nfine=701) 226 real,save :: levdat(Nfine),vmrdat(Nfine) 227 real :: vmrch4(ngrid,nlayer) ! vmr ch4 from vmrch4_proffix 228 229 REAL dtlw_nu(nlayer,L_NSPECTI) ! heating rate (K/s) due to LW in spectral bands 230 REAL dtsw_nu(nlayer,L_NSPECTV) ! heating rate (K/s) due to SW in spectral bands 231 213 232 ! local variable 233 REAL dpp ! intermediate 234 214 235 integer ok ! status (returned by NetCDF functions) 215 236 216 integer igcm_generic_ vap, igcm_generic_ice! index of the vap and ice of generic_tracer217 logical call_ice_ vap_generic ! to call only one time the ice/vap pair of a tracer237 integer igcm_generic_gas, igcm_generic_ice! index of the vap and ice of generic_tracer 238 logical call_ice_gas_generic ! to call only one time the ice/vap pair of a tracer 218 239 real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic 219 240 !$OMP THREADPRIVATE(metallicity) … … 327 348 ENDIF 328 349 329 #ifndef MESOSCALE330 350 if (is_master) call system('rm -f surf_vals_long.out') 331 #endif332 351 333 352 call su_aer_radii(ngrid,nlayer,reffrad,nueffrad) … … 434 453 435 454 436 if(varfixed .and. generic_condensation)then455 if(varfixed .and. (generic_condensation .or. methane .or. carbox))then 437 456 write(*,*) "Deep generic tracer vapor mixing ratio ? (no effect if negative) " 438 457 qvap_deep=-1. ! default value … … 504 523 endif 505 524 525 526 !----------------------------------------------------------------------- 527 ! Prepare CH4 mixing ratio for radiative transfer 528 IF (methane) then 529 vmrch4(:,:)=0. 530 531 if (ch4fix) then 532 if (vmrch4_proffix) then 533 !! Interpolate on the model vertical grid 534 do ig=1,ngrid 535 CALL interp_line(levdat,vmrdat,Nfine, & 536 zzlay(ig,:)/1000.,vmrch4(ig,:),nlayer) 537 enddo 538 else 539 vmrch4(:,:)=vmrch4fix 540 endif 541 else 542 vmrch4(:,:)=pq(:,:,igcm_ch4_gas)*100.* & 543 mmol(igcm_n2)/mmol(igcm_ch4_gas) 544 endif 545 ENDIF 546 547 ! Prepare NON LTE correction in Pluto atmosphere 548 IF (nlte) then 549 CALL nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,vmrch4,& 550 eps_nlte_sw23,eps_nlte_sw33,eps_nlte_lw) 551 ENDIF 552 ! Net atmospheric radiative cooling rate from C2H2 (K.s-1): 553 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 554 ! dtlw_hcn_c2h2=0. 555 if (cooling) then 556 ! CALL cooling_hcn_c2h2(ngrid,nlayer,pplay,& 557 ! pt,dtlw_hcn_c2h2) 558 endif 559 560 506 561 !----------------------------------------------------------------------- 507 562 do ig=1,ngrid ! Starting Big Loop over every GCM column … … 671 726 !----------------------------------------------------------------------- 672 727 673 if (generic_condensation ) then728 if (generic_condensation .or. methane .or. carbox) then 674 729 675 730 ! For now, only one GCS tracer can be both variable and radiatively active … … 679 734 do iq=1,nq 680 735 681 call generic_tracer_index(nq,iq,igcm_generic_ vap,igcm_generic_ice,call_ice_vap_generic)682 683 if (call_ice_ vap_generic) then ! to call only one time the ice/vap pair of a tracer736 call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic) 737 738 if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer 684 739 685 740 if(varactive)then 686 741 687 i_var=igcm_generic_ vap742 i_var=igcm_generic_gas 688 743 do l=1,nlayer 689 744 qvar(2*l) = pq(ig,nlayer+1-l,i_var) … … 731 786 !----------------------------------------------------------------------- 732 787 733 if (.not. generic_condensation) then788 if (.not. (generic_condensation .or. methane .or. carbox)) then 734 789 do k=1,L_LEVELS 735 790 qvar(k) = 1.0D-7 … … 741 796 ! IMPORTANT: Now convert from kg/kg to mol/mol. 742 797 do k=1,L_LEVELS 743 if (generic_condensation ) then798 if (generic_condensation .or. methane .or. carbox) then 744 799 do iq=1,nq 745 call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic) 746 if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer 747 748 epsi_generic=constants_epsi_generic(iq) 749 750 qvar(k) = qvar(k)/(epsi_generic+qvar(k)*(1.-epsi_generic)) 800 call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic) 801 if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer 802 if(.not. varactive .or. i_var.eq.iq)then 803 804 epsi_generic=constants_epsi_generic(iq) 805 806 qvar(k) = qvar(k)/(epsi_generic+qvar(k)*(1.-epsi_generic)) 807 endif 751 808 752 809 endif … … 762 819 DO l=1,nlayer 763 820 muvarrad(2*l) = muvar(ig,nlayer+2-l) 764 muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2821 muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2 765 822 END DO 766 823 … … 965 1022 tmid,pmid,taugsurf,qvar,muvarrad) 966 1023 967 call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &968 acosz,stel_fract, 969 nflux topv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu,&970 fmnetv,f luxupv,fluxdnv,fzerov,taugsurf)1024 call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, & 1025 acosz,stel_fract,nfluxtopv,fluxtopvdn,nfluxoutv_nu,& 1026 nfluxgndv_nu,nfluxtopv_nu, & 1027 fmnetv,fmnetv_nu,fluxupv,fluxdnv,fzerov,taugsurf) 971 1028 972 1029 else ! During the night, fluxes = 0. … … 1010 1067 call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, & 1011 1068 wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu, & 1012 fmneti,f luxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)1069 fmneti,fmneti_nu,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi) 1013 1070 1014 1071 !----------------------------------------------------------------------- … … 1050 1107 1051 1108 ! Finally, the heating rates 1052 1053 1109 DO l=2,L_NLAYRAD 1054 dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) & 1055 *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1))) 1056 dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) & 1057 *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1))) 1110 ! dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) & 1111 ! *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1))) 1112 dpp = glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1))) 1113 do nw=1,L_NSPECTV 1114 dtsw_nu(L_NLAYRAD+1-l,nw)= & 1115 (fmnetv_nu(l,nw)-fmnetv_nu(l-1,nw))*dpp 1116 end do 1117 1118 ! dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) & 1119 ! *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1))) 1120 do nw=1,L_NSPECTI 1121 dtlw_nu(L_NLAYRAD+1-l,nw)= & 1122 (fmneti_nu(l,nw)-fmneti_nu(l-1,nw))*dpp 1123 end do 1058 1124 END DO 1059 1125 1060 1126 ! These are values at top of atmosphere 1061 dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) & 1062 *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2))) 1063 dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) & 1064 *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2))) 1127 ! dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) & 1128 ! *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2))) 1129 ! dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) & 1130 ! *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2))) 1131 dpp = g/(cpp*scalep*(plevrad(3)-plevrad(1))) 1132 do nw=1,L_NSPECTV 1133 dtsw_nu(L_NLAYRAD,nw)= & 1134 (fmnetv_nu(1,nw)-nfluxtopv_nu(nw))*dpp 1135 end do 1136 do nw=1,L_NSPECTI 1137 dtlw_nu(L_NLAYRAD,nw)= & 1138 (fmneti_nu(1,nw)-nfluxtopi_nu(nw))*dpp 1139 end do 1065 1140 1066 1141 ! Optical thickness diagnostics (added by JVO) … … 1084 1159 endif 1085 1160 1161 ! ********************************************************** 1162 ! NON NLTE correction in Pluto atmosphere 1163 ! And conversion of LW spectral heating rates to total rates 1164 ! ********************************************************** 1165 1166 if (.not.nlte) then 1167 eps_nlte_sw23(ig,:) =1. ! IF no NLTE 1168 eps_nlte_sw33(ig,:) =1. ! IF no NLTE 1169 eps_nlte_lw(ig,:) =1. ! IF no NLTE 1170 endif 1171 1172 do l=1,nlayer 1173 1174 !LW 1175 dtlw(ig,l) =0. 1176 ! dtlw_co(ig,l) =0. ! only for diagnostic 1177 do nw=1,L_NSPECTI 1178 ! wewei : wavelength in micrometer 1179 if ((wavei(nw).gt.6.).and.(wavei(nw).lt.9)) then 1180 dtlw_nu(l,nw)=dtlw_nu(l,nw)*eps_nlte_lw(ig,l) 1181 else 1182 !dtlw_nu(l,nw)=1.*dtlw_nu(l,nw) ! no CO correction (Strobbel 1996) 1183 dtlw_nu(l,nw)=0.33*dtlw_nu(l,nw) ! CO correction (Strobbel 1996) 1184 ! dtlw_co(ig,l)=dtlw_co(ig,l)+ dtlw_nu(l,nw) ! diagnostic 1185 end if 1186 dtlw(ig,l)=dtlw(ig,l)+ dtlw_nu(l,nw) !average now on each wavelength 1187 end do 1188 ! adding c2h2 if cooling active 1189 ! dtlw(ig,l)=dtlw(ig,l)+dtlw_hcn_c2h2(ig,l) 1190 1191 !SW 1192 dtsw(ig,l) =0. 1193 1194 if (strobel) then 1195 1196 do nw=1,L_NSPECTV 1197 if ((wavev(nw).gt.2).and.(wavev(nw).lt.2.6)) then 1198 dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw23(ig,l) 1199 elseif ((wavev(nw).gt.3).and.(wavev(nw).lt.3.6)) then 1200 dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw33(ig,l) 1201 else 1202 dtsw_nu(l,nw)=dtsw_nu(l,nw) 1203 end if 1204 dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw) 1205 end do 1206 1207 else ! total heating rates multiplied by nlte coef 1208 1209 do nw=1,L_NSPECTV 1210 dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw23(ig,l) 1211 dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw) 1212 enddo 1213 1214 endif 1215 1216 1217 end do 1218 ! ********************************************************** 1219 1086 1220 1087 1221 !----------------------------------------------------------------------- 1088 1222 end do ! End of big loop over every GCM column. 1089 1223 !----------------------------------------------------------------------- 1090 1091 1224 1092 1225 -
trunk/LMDZ.PLUTO/libf/phypluto/condensation_generic_mod.F90
r3184 r3275 18 18 ! ------- 19 19 ! Calculates large-scale condensation of generic tracer "tname". 20 ! By convention, tname ends with the suffix "_ vap", as it represents the20 ! By convention, tname ends with the suffix "_gas", as it represents the 21 21 ! gas phase of the generic tracer 22 22 ! … … 39 39 REAL, intent(in) :: pdt(ngrid,nlayer) ! physical temperature tendency (K/s) 40 40 REAL, intent(in) :: pdq(ngrid,nlayer,nq) ! physical tracer tendency (K/s) 41 ! CHARACTER(*), intent(in) :: tname_ vap ! name of the tracer we consider. BY convention, it ends with _vap!!!41 ! CHARACTER(*), intent(in) :: tname_gas ! name of the tracer we consider. BY convention, it ends with _gas !!! 42 42 REAL, intent(out) :: pdtlsc(ngrid,nlayer) ! incrementation de la temperature (K) 43 43 REAL, intent(out) :: pdqvaplsc(ngrid,nlayer,nq) ! incrementation de la vapeur du traceur … … 49 49 REAL, SAVE :: qvap_deep ! deep mixing ratio of vapor when simulating bottom less planets 50 50 REAL, SAVE :: qvap_top ! top mixing ratio of vapor when simulating bottom less planets 51 logical, save :: perfect_ vap_profile52 !$OMP THREADPRIVATE(metallicity, qvap_deep, qvap_top, perfect_ vap_profile)51 logical, save :: perfect_gas_profile 52 !$OMP THREADPRIVATE(metallicity, qvap_deep, qvap_top, perfect_gas_profile) 53 53 54 54 ! Local variables 55 55 56 56 ! to call only one time the ice/vap pair of a tracer 57 logical call_ice_ vap_generic57 logical call_ice_gas_generic 58 58 59 59 INTEGER i, k , nn, iq … … 69 69 real zt(ngrid),local_p,psat_tmp,dlnpsat_tmp,Lcp,zqs_temp,zdqs 70 70 real zqs_temp_1, zqs_temp_2, zqs_temp_3 71 integer igcm_generic_ vap, igcm_generic_ice! index of the vap and ice of generic_tracer71 integer igcm_generic_gas, igcm_generic_ice! index of the vap and ice of generic_tracer 72 72 ! CHARACTER(len=*) :: tname_ice 73 73 ! evaporation calculations … … 96 96 write(*,*) " qvap_top = ",qvap_top 97 97 98 write(*,*) " perfect_ vap_profile ? "99 perfect_ vap_profile=.false. ! default value100 call getin_p("perfect_ vap_profile",perfect_vap_profile)101 write(*,*) " perfect_ vap_profile = ",perfect_vap_profile98 write(*,*) " perfect_gas_profile ? " 99 perfect_gas_profile=.false. ! default value 100 call getin_p("perfect_gas_profile",perfect_gas_profile) 101 write(*,*) " perfect_gas_profile = ",perfect_gas_profile 102 102 103 103 firstcall = .false. … … 115 115 do iq=1,nq 116 116 117 call generic_tracer_index(nq,iq,igcm_generic_ vap,igcm_generic_ice,call_ice_vap_generic)118 119 if(call_ice_ vap_generic) then ! to call only one time the ice/vap pair of a tracer117 call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic) 118 119 if(call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer 120 120 m=constants_mass(iq) 121 delta_ vapH=constants_delta_vapH(iq)121 delta_gasH=constants_delta_gasH(iq) 122 122 Tref=constants_Tref(iq) 123 123 Pref=constants_Pref(iq) … … 139 139 ! zt(i)=15. ! check too low temperatures 140 140 endif 141 zx_q(i) = pq(i,k,igcm_generic_ vap)+pdq(i,k,igcm_generic_vap)*ptimestep141 zx_q(i) = pq(i,k,igcm_generic_gas)+pdq(i,k,igcm_generic_gas)*ptimestep 142 142 143 143 call Psat_generic(zt(i),local_p,metallicity,psat_tmp,zqs_temp) … … 189 189 190 190 !Tendances de t et q 191 pdqvaplsc(1:ngrid,k,igcm_generic_ vap) = - zcond(1:ngrid)192 pdqliqlsc(1:ngrid,k,igcm_generic_ice) = - pdqvaplsc(1:ngrid,k,igcm_generic_ vap)191 pdqvaplsc(1:ngrid,k,igcm_generic_gas) = - zcond(1:ngrid) 192 pdqliqlsc(1:ngrid,k,igcm_generic_ice) = - pdqvaplsc(1:ngrid,k,igcm_generic_gas) 193 193 pdtlsc(1:ngrid,k) = pdtlsc(1:ngrid,k) + pdqliqlsc(1:ngrid,k,igcm_generic_ice)*Lcp 194 194 195 195 Enddo ! k= nlayer, 1, -1 196 196 197 if ((perfect_ vap_profile) .and. (ngrid.eq.1)) then198 ! perfect_ vap_profile is a mode that should a priori only be used in 1D:197 if ((perfect_gas_profile) .and. (ngrid.eq.1)) then 198 ! perfect_gas_profile is a mode that should a priori only be used in 1D: 199 199 ! as it is written below, it aims to force the vap profile: 200 200 ! - below condensation, it is forced to qvap_deep … … 202 202 ! - above the cold trap, it is forced to the value allowed by the cold trap 203 203 204 ! perfect_ vap_profile can be customed as you want204 ! perfect_gas_profile can be customed as you want 205 205 206 206 tau = 10. * ptimestep ! tau must not be lower than the physical step time. … … 223 223 if (k_cold_trap .lt. nlayer) then 224 224 do k=k_cold_trap+1,nlayer 225 pdqvaplsc(1,k,igcm_generic_ vap) = (pq(1,k_cold_trap,igcm_generic_vap) - pq(1,k,igcm_generic_vap))/tau - pdq(1,k,igcm_generic_vap)225 pdqvaplsc(1,k,igcm_generic_gas) = (pq(1,k_cold_trap,igcm_generic_gas) - pq(1,k,igcm_generic_gas))/tau - pdq(1,k,igcm_generic_gas) 226 226 enddo 227 227 endif … … 231 231 call Psat_generic(zt(1),pplay(1,k),metallicity,psat_tmp,zqs_temp) 232 232 if (zqs_temp .gt. qvap_deep) then 233 pdqvaplsc(1,k,igcm_generic_ vap) = (qvap_deep - pq(1,k,igcm_generic_vap))/tau - pdq(1,k,igcm_generic_vap)233 pdqvaplsc(1,k,igcm_generic_gas) = (qvap_deep - pq(1,k,igcm_generic_gas))/tau - pdq(1,k,igcm_generic_gas) 234 234 endif 235 235 if (zqs_temp .lt. qvap_deep) then 236 pdqvaplsc(1,k,igcm_generic_ vap) = (0.99*zqs_temp - pq(1,k,igcm_generic_vap))/tau - pdq(1,k,igcm_generic_vap)236 pdqvaplsc(1,k,igcm_generic_gas) = (0.99*zqs_temp - pq(1,k,igcm_generic_gas))/tau - pdq(1,k,igcm_generic_gas) 237 237 endif 238 238 enddo … … 246 246 ! brings lower generic vapor ratio to a fixed value. 247 247 ! tau is in seconds and must not be lower than the physical step time. 248 pdqvaplsc(1:ngrid,1,igcm_generic_ vap) = (qvap_deep - pq(1:ngrid,1,igcm_generic_vap))/tau - pdq(1:ngrid,1,igcm_generic_vap)248 pdqvaplsc(1:ngrid,1,igcm_generic_gas) = (qvap_deep - pq(1:ngrid,1,igcm_generic_gas))/tau - pdq(1:ngrid,1,igcm_generic_gas) 249 249 endif 250 250 if (qvap_top >= 0.) then … … 252 252 ! brings lower generic vapor ratio to a fixed value. 253 253 ! tau is in seconds and must not be lower than the physical step time. 254 pdqvaplsc(1:ngrid,nlayer,igcm_generic_ vap) = (qvap_top - pq(1:ngrid,nlayer,igcm_generic_vap))/tau - pdq(1:ngrid,nlayer,igcm_generic_vap)254 pdqvaplsc(1:ngrid,nlayer,igcm_generic_gas) = (qvap_top - pq(1:ngrid,nlayer,igcm_generic_gas))/tau - pdq(1:ngrid,nlayer,igcm_generic_gas) 255 255 endif 256 endif !if(call_ice_ vap_generic)256 endif !if(call_ice_gas_generic) 257 257 enddo ! iq=1,nq 258 258 -
trunk/LMDZ.PLUTO/libf/phypluto/convadj.F
r3195 r3275 304 304 end do 305 305 ! do l = 1, nlay 306 ! print*,'zq(ig,:,vap) = ',zq(i,l,igcm_h2o_ vap)306 ! print*,'zq(ig,:,vap) = ',zq(i,l,igcm_h2o_gas) 307 307 ! end do 308 308 ! do l = 1, nlay 309 ! print*,'zq2(ig,:,vap) = ',zq2(i,l,igcm_h2o_ vap)309 ! print*,'zq2(ig,:,vap) = ',zq2(i,l,igcm_h2o_gas) 310 310 ! end do 311 ! print*,'zqm(vap) = ',zqm(igcm_h2o_ vap)311 ! print*,'zqm(vap) = ',zqm(igcm_h2o_gas) 312 312 print*,'jadrs=',jadrs 313 313 -
trunk/LMDZ.PLUTO/libf/phypluto/dyn1d
-
Property
svn:ignore
set to
*.mod
-
Property
svn:ignore
set to
-
trunk/LMDZ.PLUTO/libf/phypluto/generic_cloud_common_h.F90
r3184 r3275 16 16 implicit none 17 17 real, save :: m ! molecular mass of the specie (g/mol) 18 real, save :: delta_ vapH ! Enthalpy of vaporization (J/mol)18 real, save :: delta_gasH ! Enthalpy of vaporization (J/mol) 19 19 real,save :: Tref ! Ref temperature for Clausis-Clapeyron (K) 20 20 real,save :: Pref ! Reference pressure for Clausius Clapeyron (Pa) … … 46 46 if (trim(specname) .eq. "Mg") then 47 47 print*,"Loading data for Mg" 48 delta_ vapH = 521.7E3 !J/mol48 delta_gasH = 521.7E3 !J/mol 49 49 Tref = 2303.0 !K 50 50 Pref = 1.0E5 !in Pa 51 51 m = 140.6931 52 RLVTT_generic = delta_ vapH/(m/1000.)52 RLVTT_generic = delta_gasH/(m/1000.) 53 53 metallicity_coeff=1.0*log(10.0) 54 54 55 55 else if (trim(specname) .eq. "Na") then 56 56 print*,"Loading data for Na" 57 delta_ vapH = 265.9E357 delta_gasH = 265.9E3 58 58 Tref = 1624.0 59 59 Pref = 1.0E5 60 60 m = 78.04454 !Na2S 61 RLVTT_generic = delta_ vapH/(m/1000.)61 RLVTT_generic = delta_gasH/(m/1000.) 62 62 metallicity_coeff=0.5*log(10.0) 63 63 64 64 else if (trim(specname) .eq. "Fe") then 65 65 print*,"Loading data for Fe" 66 delta_ vapH = 401.9E366 delta_gasH = 401.9E3 67 67 Tref = 2903.0 68 68 Pref = 1.0E5 69 69 m = 55.8450 70 RLVTT_generic = delta_ vapH/(m/1000.)70 RLVTT_generic = delta_gasH/(m/1000.) 71 71 metallicity_coeff=0.0*log(10.0) 72 72 73 73 else if (trim(specname) .eq. "Cr") then 74 74 print*,"Loading data for Cr" 75 delta_ vapH = 394.2E375 delta_gasH = 394.2E3 76 76 Tref = 2749.0 77 77 Pref = 1.0E5 78 78 m = 51.9961 79 RLVTT_generic = delta_ vapH/(m/1000.)79 RLVTT_generic = delta_gasH/(m/1000.) 80 80 metallicity_coeff=0.0*log(10.0) 81 81 82 82 else if (trim(specname) .eq. "KCl") then 83 83 print*,"Loading data for KCl" 84 delta_ vapH = 217.9E384 delta_gasH = 217.9E3 85 85 Tref = 1495.0 86 86 Pref = 1.0E5 87 87 metallicity_coeff=0.0*log(10.0) 88 88 m = 74.5498 89 RLVTT_generic = delta_ vapH/(m/1000.)89 RLVTT_generic = delta_gasH/(m/1000.) 90 90 else if (trim(specname) .eq. "Mn") then 91 91 print*,"Loading data for Mn" 92 delta_ vapH = 455.8E392 delta_gasH = 455.8E3 93 93 Tref = 2064.0 94 94 Pref = 1.0E5 95 95 metallicity_coeff=1.0*log(10.0) 96 96 m = 87.003049 97 RLVTT_generic = delta_ vapH/(m/1000.)97 RLVTT_generic = delta_gasH/(m/1000.) 98 98 else if (trim(specname) .eq. "Zn") then 99 99 print*,"Loading data for Zn" 100 delta_ vapH = 303.9E3100 delta_gasH = 303.9E3 101 101 Tref = 1238.0 102 102 Pref = 1.0E5 103 103 metallicity_coeff=1.0*log(10.0) 104 104 m = 97.445 105 RLVTT_generic = delta_ vapH/(m/1000.)105 RLVTT_generic = delta_gasH/(m/1000.) 106 106 else 107 107 print*,"Unknow species (not in Mg, Fe, Na, KCl, Cr, Mn or Zn)" … … 138 138 if (index(table_traceurs_line,specname) /= 0) then 139 139 140 if (index(table_traceurs_line,'delta_ vapH=' ) /= 0) then141 read(table_traceurs_line(index(table_traceurs_line,'delta_ vapH=')+len('delta_vapH='):),*) delta_vapH142 print*, 'delta_ vapH ', delta_vapH140 if (index(table_traceurs_line,'delta_gasH=' ) /= 0) then 141 read(table_traceurs_line(index(table_traceurs_line,'delta_gasH=')+len('delta_gasH='):),*) delta_gasH 142 print*, 'delta_gasH ', delta_gasH 143 143 end if 144 144 if (index(table_traceurs_line,'Tref=' ) /= 0) then … … 178 178 close(117) 179 179 180 RLVTT_generic=delta_ vapH/(m/1000.)180 RLVTT_generic=delta_gasH/(m/1000.) 181 181 182 182 write(*,*) 'RLVTT_generic', RLVTT_generic … … 190 190 !============================================================================ 191 191 ! Clausius-Clapeyron relation : 192 ! d(Ln(psat))/dT = delta_ vapH/(RT^2)193 ! ->Psat = Pref * exp(-delta_ vapH/R*(1/T-1/Tref))192 ! d(Ln(psat))/dT = delta_gasH/(RT^2) 193 ! ->Psat = Pref * exp(-delta_gasH/R*(1/T-1/Tref)) 194 194 ! 195 195 ! Authors … … 203 203 real, intent(out) :: psat,qsat !saturation vapor pressure (Pa) and mass mixing ratio at saturation (kg/kg) of the layer 204 204 205 psat = pref*exp(-delta_ vapH/(r*mugaz/1000.)*(1/T-1/Tref)-metallicity_coeff*metallicity) ! in Pa (because pref in Pa)205 psat = pref*exp(-delta_gasH/(r*mugaz/1000.)*(1/T-1/Tref)-metallicity_coeff*metallicity) ! in Pa (because pref in Pa) 206 206 207 207 if (psat .gt. p) then … … 217 217 !=============================================================================== 218 218 ! Compute dqsat = L/cp* d(qsat)/dT and d(Ln(psat)) = L/cp*d(Ln(psat))/dT 219 ! we have d(ln(psat))/dT = delta_ vapH/R*(1/T^2) for clausius-clapeyron219 ! we have d(ln(psat))/dT = delta_gasH/R*(1/T^2) for clausius-clapeyron 220 220 ! r*mugaz/1000. is the perfect gaz constant in the computation of "dummy" 221 221 ! Authors … … 235 235 real :: dummy ! used to store d(ln(psat))/dT 236 236 237 dummy = delta_ vapH/((r*mugaz/1000.)*(T**2))237 dummy = delta_gasH/((r*mugaz/1000.)*(T**2)) 238 238 239 239 if (psat .gt. p) then … … 270 270 271 271 272 Tsat = 1/(1/Tref - (r*mugaz/1000.)/delta_ vapH*Log(p/Pref))272 Tsat = 1/(1/Tref - (r*mugaz/1000.)/delta_gasH*Log(p/Pref)) 273 273 274 274 end subroutine Tsat_generic -
trunk/LMDZ.PLUTO/libf/phypluto/generic_tracer_index_mod.F90
r3184 r3275 6 6 ! for both the vap & ice parts of the tracer 7 7 ! 8 ! call_ice_ vap_generic is a boolean which tells if the tracer is condensable9 ! and allows, if it is condensable, to call the vap/ice pair 8 ! call_ice_gas_generic is a boolean which tells if the tracer is condensable 9 ! and allows, if it is condensable, to call the vap/ice pair 10 10 ! 11 11 ! Noé Clément & Lucas Teinturier (2022) … … 17 17 contains 18 18 19 subroutine generic_tracer_index(nq,iq,igcm_generic_ vap,igcm_generic_ice,call_ice_vap_generic)19 subroutine generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic) 20 20 21 21 USE tracer_h, only: noms, is_condensable 22 integer nq,iq, igcm_generic_ vap, igcm_generic_ice, i_shortname, ii23 logical call_ice_ vap_generic22 integer nq,iq, igcm_generic_gas, igcm_generic_ice, i_shortname, ii 23 logical call_ice_gas_generic 24 24 character(len=10) :: shortname 25 if((is_condensable(iq)==1) .and. (index(noms(iq)," vap") .ne. 0) &25 if((is_condensable(iq)==1) .and. (index(noms(iq),"gas") .ne. 0) & 26 26 .and. (index(noms(iq),"h2o") .eq. 0) .and. (index(noms(iq),"n2") .eq. 0)) then 27 ! Let's get the index of our tracers (we look for igcm _generic_ vapand igcm_generic_ice)27 ! Let's get the index of our tracers (we look for igcm _generic_gas and igcm_generic_ice) 28 28 29 igcm_generic_ vap= iq29 igcm_generic_gas = iq 30 30 igcm_generic_ice = -1 31 31 … … 46 46 ! ! endif 47 47 48 ! ! call_ice_ vap_generic = .true.48 ! ! call_ice_gas_generic = .true. 49 49 i_shortname = index(noms(iq),'_')-1 50 50 shortname=noms(iq)(1:i_shortname) … … 59 59 endif 60 60 ! if we didn't find it before, we look after the vap tracer in pq 61 if (igcm_generic_ice .eq. -1 .and. (iq<nq)) then 61 if (igcm_generic_ice .eq. -1 .and. (iq<nq)) then 62 62 do ii=iq+1,nq 63 if (index(noms(ii),adjustl(trim(shortname))) .ne. 0) then 63 if (index(noms(ii),adjustl(trim(shortname))) .ne. 0) then 64 64 igcm_generic_ice=ii 65 65 exit … … 67 67 enddo 68 68 endif 69 call_ice_ vap_generic = .true.69 call_ice_gas_generic = .true. 70 70 else 71 71 72 call_ice_ vap_generic = .false.72 call_ice_gas_generic = .false. 73 73 74 74 end if ! if((is_condensable(iq)==1) .and. (index(noms(iq),"vap") .ne. 0)) -
trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90
r3258 r3275 883 883 " albmin_ch4 = ",albmin_ch4 884 884 885 886 if (is_master)write(*,*)trim(rname)//& 887 "fixed gaseous CH4 mixing ratio for rad transfer ?" 888 ch4fix=.false. ! default value 889 call getin_p("ch4fix",ch4fix) 890 if (is_master)write(*,*)trim(rname)//& 891 " ch4fix = ",ch4fix 892 if (is_master)write(*,*)trim(rname)//& 893 "fixed gaseous CH4 mixing ratio for rad transfer ?" 894 vmrch4fix=0.5 ! default value 895 call getin_p("vmrch4fix",vmrch4fix) 896 if (is_master)write(*,*)trim(rname)//& 897 " vmrch4fix = ",vmrch4fix 898 vmrch4_proffix=.false. ! default value 899 call getin_p("vmrch4_proffix",vmrch4_proffix) 900 if (is_master)write(*,*)trim(rname)//& 901 " vmrch4_proffix = ",vmrch4_proffix 902 903 885 904 !*********** CO ********************************* 886 905 -
trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90
r3258 r3275 127 127 IF (.NOT. allocated(is_rgcs)) allocate(is_rgcs(nq)) !LT 128 128 IF (.NOT. allocated(constants_mass)) allocate(constants_mass(nq)) 129 IF (.NOT. allocated(constants_delta_ vapH)) allocate(constants_delta_vapH(nq))129 IF (.NOT. allocated(constants_delta_gasH)) allocate(constants_delta_gasH(nq)) 130 130 IF (.NOT. allocated(constants_Tref)) allocate(constants_Tref(nq)) 131 131 IF (.NOT. allocated(constants_Pref)) allocate(constants_Pref(nq)) … … 157 157 is_rgcs(:) = 0 158 158 constants_mass(:)=0 159 constants_delta_ vapH(:)=0159 constants_delta_gasH(:)=0 160 160 constants_Tref(:)=0 161 161 constants_Pref(:)=0 … … 229 229 write(*,*) 'Tracer ',count,' = n2' 230 230 endif 231 if (noms(iq).eq."ch4_ vap") then231 if (noms(iq).eq."ch4_gas") then 232 232 igcm_ch4_gas=iq 233 233 mmol(igcm_ch4_gas)=16. … … 243 243 write(*,*) 'Tracer ',count,' = ch4 ice' 244 244 endif 245 if (noms(iq).eq."co_ vap") then245 if (noms(iq).eq."co_gas") then 246 246 igcm_co_gas=iq 247 247 mmol(igcm_co_gas)=28. … … 376 376 call specie_parameters_table(noms(iq)(1:len(trim(noms(iq)))-4)) 377 377 constants_mass(iq)=m 378 constants_delta_ vapH(iq)=delta_vapH378 constants_delta_gasH(iq)=delta_gasH 379 379 constants_Tref(iq)=Tref 380 380 constants_Pref(iq)=Pref -
trunk/LMDZ.PLUTO/libf/phypluto/mass_redistribution_mod.F90
r3184 r3275 7 7 pu,pv,pdt,pdtsrf,pdq,pdu,pdv,pdmassmr, & 8 8 pdtmr,pdtsrfmr,pdpsrfmr,pdumr,pdvmr,pdqmr,pdqsmr) 9 9 10 10 use radcommon_h, only: glat 11 11 USE planete_mod, only: bp 12 12 use comcstfi_mod, only: g 13 13 14 14 IMPLICIT NONE 15 15 !======================================================================= … … 22 22 ! 23 23 ! input: 24 ! ----- 24 ! ----- 25 25 ! ngrid nombre de points de verticales 26 26 ! (toutes les boucles de la physique sont au 27 27 ! moins vectorisees sur ngrid) 28 28 ! nlayer nombre de couches 29 ! pplay(ngrid,nlayer) Pressure levels 30 ! pplev(ngrid,nlayer+1) Pressure levels 29 ! pplay(ngrid,nlayer) Pressure levels 30 ! pplev(ngrid,nlayer+1) Pressure levels 31 31 ! nq Number of tracers 32 32 ! … … 35 35 ! pu,pv (ngrid,nlayer) wind velocity (m/s) 36 36 ! 37 ! 37 ! 38 38 ! pdX physical tendency of X before mass redistribution 39 39 ! … … 45 45 ! pdXmr(ngrid) physical tendency of X after mass redistribution 46 46 ! 47 ! 47 ! 48 48 ! 49 49 !======================================================================= … … 55 55 ! Arguments : 56 56 ! --------- 57 INTEGER,INTENT(IN) :: ngrid, nlayer, nq 57 INTEGER,INTENT(IN) :: ngrid, nlayer, nq 58 58 REAL,INTENT(IN) :: ptimestep 59 59 REAL,INTENT(IN) :: pcapcal(ngrid) … … 80 80 REAL zzu(nlayer),zzv(nlayer) 81 81 REAL zzq(nlayer,nq),zzt(nlayer) 82 ! Dummy variables 82 ! Dummy variables 83 83 INTEGER n,l,ig,iq 84 84 REAL zdtsig(ngrid,nlayer) … … 88 88 REAL zq1(nlayer) 89 89 REAL ztsrf(ngrid) 90 REAL ztm(nlayer+1) 90 REAL ztm(nlayer+1) 91 91 REAL zum(nlayer+1) , zvm(nlayer+1) 92 92 REAL zqm(nlayer+1,nq),zqm1(nlayer+1) … … 103 103 ! 104 104 IF (firstcall) THEN 105 firstcall=.false. 105 firstcall=.false. 106 106 ENDIF 107 107 ! 108 108 !====================================================================== 109 ! Calcul of h2o condensation 109 ! Calcul of h2o condensation 110 110 ! ============================================================ 111 ! 111 ! 112 112 ! Used variable : 113 113 ! pdmassmr : air Mass added to the atmosphere in each layer per unit time (kg/m2/s) … … 135 135 ! Changing pressure at the surface: 136 136 ! """""""""""""""""""""""""""""""""""" 137 137 138 zmassboil(1:ngrid) = 0 ! AF: TODO: update this for pluto 139 pdtsrfmr(1:ngrid) = 0 140 138 141 pdpsrfmr(1:ngrid) = (zdmass_sum(1:ngrid,1)+zmassboil(1:ngrid))*g 139 142 … … 151 154 152 155 ! *************************************************************** 153 ! Correction to account for redistribution between sigma or hybrid 156 ! Correction to account for redistribution between sigma or hybrid 154 157 ! layers when changing surface pressure 155 158 ! zzx quantities have dimension (nlayer) to speed up calculation … … 174 177 ! Ehouarn: but for now leave things as before 175 178 zmflux(l+1) = zmflux(l) + pdmassmr(ig,l) - (bp(l)-bp(l+1))*(zdmass_sum(ig,1)+zmflux(1)) 176 ! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld 179 ! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld 177 180 if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0. 178 181 END DO 179 182 180 183 ! Mass of each layer 181 ! ------------------ 184 ! ------------------ 182 185 zzmass(1:nlayer)=zmass(ig,1:nlayer)*(1.+pdpsrfmr(ig)*ptimestep/pplev(ig,1)) 183 186 … … 186 189 ! """""""""""""""""""""""""""""""" 187 190 188 ! averaging operator for TRANSPORT 191 ! averaging operator for TRANSPORT 189 192 ! """""""""""""""""""""""""""""""" 190 193 ! Value transfert at the surface interface when condensation 191 194 ! sublimation: 192 195 ztm(1) = ztsrf(ig) 193 zum(1) = 0. 194 zvm(1) = 0. 196 zum(1) = 0. 197 zvm(1) = 0. 195 198 zqm(1,1:nq)=0. ! most tracer do not condense ! 196 !! if (water) zqm(1,igcm_h2o_ vap)=1. ! flux is 100% h2o at surface197 199 !! if (water) zqm(1,igcm_h2o_gas)=1. ! flux is 100% h2o at surface 200 198 201 ! Van Leer scheme: 199 202 w(1:nlayer+1)=-zmflux(1:nlayer+1)*ptimestep 200 call vl1d(nlayer,zzt,2.,zzmass,w,ztm) 201 call vl1d(nlayer,zzu,2.,zzmass,w,zum) 202 call vl1d(nlayer,zzv,2.,zzmass,w,zvm) 203 call vl1d(nlayer,zzt,2.,zzmass,w,ztm) 204 call vl1d(nlayer,zzu,2.,zzmass,w,zum) 205 call vl1d(nlayer,zzv,2.,zzmass,w,zvm) 203 206 do iq=1,nq 204 207 zq1(1:nlayer)=zzq(1:nlayer,iq) … … 214 217 215 218 ! Surface condensation affects low winds 216 if (zmflux(1).lt.0) then 219 if (zmflux(1).lt.0) then 217 220 zum(1)= zzu(1) * (w(1)/zzmass(1)) 218 221 zvm(1)= zzv(1) * (w(1)/zzmass(1)) … … 222 225 end if 223 226 end if 224 225 ztm(nlayer+1)= zzt(nlayer) ! should not be used, but... 227 228 ztm(nlayer+1)= zzt(nlayer) ! should not be used, but... 226 229 zum(nlayer+1)= zzu(nlayer) ! should not be used, but... 227 230 zvm(nlayer+1)= zzv(nlayer) ! should not be used, but... 228 231 zqm(nlayer+1,1:nq)= zzq(nlayer,1:nq) 229 230 ! Tendencies on T, U, V, Q 232 233 ! Tendencies on T, U, V, Q 231 234 ! """""""""""""""""""""""" 232 235 DO l=1,nlayer 233 236 234 237 ! Tendencies on T 235 238 pdtmr(ig,l) = (1/zzmass(l)) * & … … 253 256 enddo 254 257 255 END DO ! loop on ig 258 END DO ! loop on ig 256 259 257 260 CONTAINS … … 260 263 SUBROUTINE vl1d(llm,q,pente_max,zzmass,w,qm) 261 264 ! 262 ! 265 ! 263 266 ! Operateur de moyenne inter-couche pour calcul de transport type 264 267 ! Van-Leer " pseudo amont " dans la verticale … … 270 273 ! 271 274 ! -------------------------------------------------------------------- 272 275 273 276 IMPLICIT NONE 274 277 … … 280 283 REAL w(llm+1) 281 284 ! 282 ! Local 285 ! Local 283 286 ! --------- 284 287 ! … … 287 290 real dzq(llm),dzqw(llm),adzqw(llm),dzqmax 288 291 real sigw, Mtot, MQtot 289 integer m 290 ! integer ismax,ismin 291 292 293 ! On oriente tout dans le sens de la pression 292 integer m 293 ! integer ismax,ismin 294 295 296 ! On oriente tout dans le sens de la pression 294 297 ! W > 0 WHEN DOWN !!!!!!!!!!!!! 295 298 … … 347 350 stop 348 351 end if 349 else ! if(w(l+1).lt.0) 350 m = l-1 352 else ! if(w(l+1).lt.0) 353 m = l-1 351 354 Mtot = zzmass(m+1) 352 355 MQtot = zzmass(m+1)*q(m+1) … … 365 368 else 366 369 qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*qm(1)) 367 end if 370 end if 368 371 end if 369 372 enddo … … 373 376 ! if(w(1).gt.0.) then 374 377 ! qm(1)=q(1) 375 ! else 378 ! else 376 379 ! qm(1)=0. 377 380 ! end if -
trunk/LMDZ.PLUTO/libf/phypluto/phyetat0_mod.F90
r3184 r3275 257 257 endif ! of if (nq.ge.1) 258 258 259 !!call WriteField_phy("in_phyetat0_qsurf",qsurf(1:ngrid,igcm_h2o_ vap),1)259 !!call WriteField_phy("in_phyetat0_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1) 260 260 261 261 if (startphy_file) then -
trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90
r3258 r3275 300 300 real omega(ngrid,nlayer) ! omega velocity (Pa/s, >0 when downward) 301 301 302 integer i,l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_ vap, igcm_generic_ice303 logical call_ice_ vap_generic ! to call only one time the ice/vap pair of a tracer302 integer i,l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_gas, igcm_generic_ice 303 logical call_ice_gas_generic ! to call only one time the ice/vap pair of a tracer 304 304 305 305 real zls ! Solar longitude (radians). … … 638 638 qsurf_hist(:,:)=qsurf(:,:) 639 639 640 !! call WriteField_phy("post_qsurf_hist_qsurf",qsurf(1:ngrid,igcm_h2o_ vap),1)640 !! call WriteField_phy("post_qsurf_hist_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1) 641 641 642 642 ! Initialize variable for dynamical heating and zonal wind tendency diagnostic … … 681 681 endif 682 682 683 !! call WriteField_phy("post_physdem_qsurf",qsurf(1:ngrid,igcm_h2o_ vap),1)683 !! call WriteField_phy("post_physdem_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1) 684 684 #endif 685 685 if (corrk) then … … 702 702 endif 703 703 704 !! call WriteField_phy("post_corrk_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_ vap),1)704 !! call WriteField_phy("post_corrk_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1) 705 705 ! XIOS outputs 706 706 #ifdef CPP_XIOS … … 711 711 #endif 712 712 713 !! call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_ vap),1)713 !! call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1) 714 714 715 715 write(*,*) "physiq: end of firstcall" 716 716 endif ! end of 'firstcall' 717 717 718 !! call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_ vap),1)719 !! call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_ vap))718 !! call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1) 719 !! call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas)) 720 720 721 721 if (check_physics_inputs) then … … 882 882 ! Eclipse incoming sunlight !AF24: removed 883 883 884 !! call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_ vap))885 !! call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_ vap))884 !! call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas)) 885 !! call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas)) 886 886 887 887 … … 900 900 do iq=1,nq 901 901 902 call generic_tracer_index(nq,iq,igcm_generic_ vap,igcm_generic_ice,call_ice_vap_generic)903 904 if (call_ice_ vap_generic) then ! to call only one time the ice/vap pair of a tracer902 call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic) 903 904 if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer 905 905 906 906 epsi_generic=constants_epsi_generic(iq) 907 907 908 muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_ vap))909 muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_ vap))908 muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_gas)) 909 muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_gas)) 910 910 911 911 endif … … 931 931 call callcorrk(ngrid,nlayer,pq,nq,qsurf, & 932 932 albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, & 933 tsurf,fract,dist_star,aerosol,muvar, &933 zzlay,tsurf,fract,dist_star,aerosol,muvar, & 934 934 zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, & 935 935 fluxsurfabs_sw,fluxtop_lw, & … … 1039 1039 endif ! of if (callrad) 1040 1040 1041 !! call writediagfi(ngrid,"vdifc_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_ vap))1042 !! call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_ vap))1041 !! call writediagfi(ngrid,"vdifc_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas)) 1042 !! call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas)) 1043 1043 1044 1044 … … 1137 1137 ! endif 1138 1138 1139 !! call writediagfi(ngrid,"vdifc_post_zdqsdif"," "," ",2,zdqsdif(1:ngrid,igcm_h2o_ vap))1139 !! call writediagfi(ngrid,"vdifc_post_zdqsdif"," "," ",2,zdqsdif(1:ngrid,igcm_h2o_gas)) 1140 1140 1141 1141 if (tracer) then … … 1144 1144 end if ! of if (tracer) 1145 1145 1146 !! call writediagfi(ngrid,"vdifc_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_ vap))1147 !! call writediagfi(ngrid,"vdifc_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_ vap))1146 !! call writediagfi(ngrid,"vdifc_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas)) 1147 !! call writediagfi(ngrid,"vdifc_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas)) 1148 1148 1149 1149 ! test energy conservation … … 1264 1264 dqsurf(1:ngrid,igcm_n2) = dqsurf(1:ngrid,igcm_n2) + zdqsc(1:ngrid,igcm_n2) 1265 1265 1266 !! call writediagfi(ngrid,"condense_n2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_ vap))1267 !! call writediagfi(ngrid,"condense_n2_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_ vap))1266 !! call writediagfi(ngrid,"condense_n2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas)) 1267 !! call writediagfi(ngrid,"condense_n2_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas)) 1268 1268 1269 1269 ! test energy conservation … … 1458 1458 cloudfrac(:,:)=0.0 1459 1459 do iq=1,nq 1460 call generic_tracer_index(nq,iq,igcm_generic_ vap,igcm_generic_ice,call_ice_vap_generic)1461 if (call_ice_ vap_generic) then ! to call only one time the ice/vap pair of a tracer1460 call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic) 1461 if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer 1462 1462 do l = 1, nlayer 1463 1463 do ig=1,ngrid … … 1491 1491 dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq) 1492 1492 1493 !! call writediagfi(ngrid,"callsedim_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_ vap))1493 !! call writediagfi(ngrid,"callsedim_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas)) 1494 1494 1495 1495 ! ! Test water conservation !AF24: removed … … 1497 1497 end if ! end of 'sedimentation' 1498 1498 1499 !! call writediagfi(ngrid,"mass_redist_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_ vap))1500 !! call writediagfi(ngrid,"mass_redist_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_ vap))1499 !! call writediagfi(ngrid,"mass_redist_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas)) 1500 !! call writediagfi(ngrid,"mass_redist_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas)) 1501 1501 1502 1502 ! --------------- … … 1509 1509 zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * 0 1510 1510 ! ( zdqevap(1:ngrid,1:nlayer) & 1511 ! ! + zdqrain(1:ngrid,1:nlayer,igcm_h2o_ vap) &1512 ! ! + dqmoist(1:ngrid,1:nlayer,igcm_h2o_ vap) &1511 ! ! + zdqrain(1:ngrid,1:nlayer,igcm_h2o_gas) & 1512 ! ! + dqmoist(1:ngrid,1:nlayer,igcm_h2o_gas) & 1513 1513 ! + dqvaplscale(1:ngrid,1:nlayer) ) 1514 1514 … … 1536 1536 endif 1537 1537 1538 ! call writediagfi(ngrid,"mass_redistribution_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_ vap))1539 1540 !! call writediagfi(ngrid,"slab_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_ vap))1541 !! call writediagfi(ngrid,"slab_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_ vap))1538 ! call writediagfi(ngrid,"mass_redistribution_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas)) 1539 1540 !! call writediagfi(ngrid,"slab_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas)) 1541 !! call writediagfi(ngrid,"slab_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas)) 1542 1542 1543 1543 … … 1794 1794 do iq=1,nq 1795 1795 1796 call generic_tracer_index(nq,iq,igcm_generic_ vap,igcm_generic_ice,call_ice_vap_generic)1797 1798 if (call_ice_ vap_generic) then ! to call only one time the ice/vap pair of a tracer1796 call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic) 1797 1798 if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer 1799 1799 1800 1800 do l = 1, nlayer 1801 1801 do ig=1,ngrid 1802 1802 call Psat_generic(zt(ig,l),pplay(ig,l),metallicity,psat_tmp_generic,qsat_generic(ig,l,iq)) 1803 RH_generic(ig,l,iq) = zq(ig,l,igcm_generic_ vap) / qsat_generic(ig,l,iq)1803 RH_generic(ig,l,iq) = zq(ig,l,igcm_generic_gas) / qsat_generic(ig,l,iq) 1804 1804 enddo 1805 1805 enddo … … 2364 2364 2365 2365 do iq=1,nq 2366 call generic_tracer_index(nq,iq,igcm_generic_ vap,igcm_generic_ice,call_ice_vap_generic)2367 2368 if (call_ice_ vap_generic) then ! to call only one time the ice/vap pair of a tracer2366 call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic) 2367 2368 if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer 2369 2369 2370 2370 reffrad_generic_zeros_for_wrf(:,:) = 1. … … 2373 2373 comm_TOTCLOUDFRAC(1:ngrid) = totcloudfrac(1:ngrid) !?????? 2374 2374 comm_SURFRAIN(1:ngrid) = zdqsrain_generic(1:ngrid,iq) 2375 comm_DQVAP(1:ngrid,1:nlayer) = pdq(1:ngrid,1:nlayer,igcm_generic_ vap)2375 comm_DQVAP(1:ngrid,1:nlayer) = pdq(1:ngrid,1:nlayer,igcm_generic_gas) 2376 2376 comm_DQICE(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_generic_ice) 2377 2377 ! comm_H2OICE_REFF(1:ngrid,1:nlayer) = reffrad_generic_zeros_for_wrf(1:ngrid,1:nlayer) ! for now zeros ! … … 2417 2417 ! do ig=1,ngrid 2418 2418 ! if ((tracer).and.(water)) then 2419 ! pdq(ig,:,igcm_h2o_ vap) = pdq(ig,:,igcm_h2o_vap) + lsf_dq(:)2419 ! pdq(ig,:,igcm_h2o_gas) = pdq(ig,:,igcm_h2o_gas) + lsf_dq(:) 2420 2420 ! endif 2421 2421 ! pdt(ig,:) = pdt(ig,:) + lsf_dt(:) -
trunk/LMDZ.PLUTO/libf/phypluto/sfluxi.F
r3184 r3275 1 1 module sfluxi_mod 2 2 3 3 implicit none 4 4 5 5 contains 6 6 7 7 SUBROUTINE SFLUXI(PLEV,TLEV,DTAUI,TAUCUMI,UBARI,RSFI,WNOI,DWNI, 8 8 * COSBI,WBARI,NFLUXTOPI,NFLUXTOPI_nu, 9 * FMNETI,f luxupi,fluxdni,fluxupi_nu,9 * FMNETI,fmneti_nu,fluxupi,fluxdni,fluxupi_nu, 10 10 * FZEROI,TAUGSURF) 11 11 12 12 use radinc_h, only: NTfac, NTstart, L_LEVELS, L_NSPECTI, L_NGAUSS 13 13 use radinc_h, only: L_NLAYRAD, L_NLEVRAD … … 15 15 use comcstfi_mod, only: pi 16 16 use gfluxi_mod, only: gfluxi 17 17 18 18 implicit none 19 19 20 20 integer NLEVRAD, L, NW, NG, NTS, NTT 21 21 22 22 real*8 TLEV(L_LEVELS), PLEV(L_LEVELS) 23 23 real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS) 24 24 real*8 FMNETI(L_NLAYRAD) 25 real*8 FMNETI_NU(L_NLAYRAD,L_NSPECTI) 25 26 real*8 WNOI(L_NSPECTI), DWNI(L_NSPECTI) 26 27 real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS) … … 32 33 real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) 33 34 real*8 FTOPUP 34 35 35 36 real*8 UBARI, RSFI, TSURF, BSURF, TTOP, BTOP, TAUTOP 36 37 real*8 PLANCK, PLTOP … … 38 39 real*8 FZEROI(L_NSPECTI) 39 40 real*8 taugsurf(L_NSPECTI,L_NGAUSS-1), fzero 40 41 41 42 real*8 fup_tmp(L_NSPECTI),fdn_tmp(L_NSPECTI) 42 43 real*8 PLANCKSUM,PLANCKREF 43 44 44 45 ! AB : variables for interpolation 45 46 REAL*8 C1 46 47 REAL*8 C2 47 48 REAL*8 P1 48 49 49 50 !======================================================================C 50 51 51 52 NLEVRAD = L_NLEVRAD 52 53 53 54 ! ZERO THE NET FLUXES 54 55 NFLUXTOPI = 0.0D0 55 56 56 57 DO NW=1,L_NSPECTI 57 58 NFLUXTOPI_nu(NW) = 0.0D0 58 59 DO L=1,L_NLAYRAD 59 60 FLUXUPI_nu(L,NW) = 0.0D0 61 FMNETI_nu(L,NW) = 0.0 60 62 fup_tmp(nw)=0.0D0 61 63 fdn_tmp(nw)=0.0D0 62 64 END DO 63 65 END DO 64 66 65 67 DO L=1,L_NLAYRAD 66 68 FMNETI(L) = 0.0D0 … … 68 70 FLUXDNI(L) = 0.0D0 69 71 END DO 70 72 71 73 ! WE NOW ENTER A MAJOR LOOP OVER SPECTRAL INTERVALS IN THE INFRARED 72 74 ! TO CALCULATE THE NET FLUX IN EACH SPECTRAL INTERVAL 73 75 74 76 TTOP = TLEV(2) ! JL12 why not (1) ??? 75 77 TSURF = TLEV(L_LEVELS) … … 86 88 PLANCKREF = TSURF * TSURF 87 89 PLANCKREF = sigma * PLANCKREF * PLANCKREF 88 90 89 91 DO NW=1,L_NSPECTI 90 92 ! AB : PLANCKIR(NW,NTS) is replaced by P1, the linear interpolation result for a temperature TSURF … … 93 95 PLANCKSUM = PLANCKSUM + P1 * DWNI(NW) 94 96 ENDDO 95 97 96 98 PLANCKSUM = PLANCKREF / (PLANCKSUM * Pi) 97 99 !JL12 98 100 99 101 DO 501 NW=1,L_NSPECTI 100 102 ! SURFACE EMISSIONS - INDEPENDENT OF GAUSS POINTS … … 106 108 BSURF = (1. - RSFI) * P1 * PLANCKSUM 107 109 PLTOP = (1.0D0 - C2) * PLANCKIR(NW,NTT) + C2*PLANCKIR(NW,NTT+1) 108 110 109 111 ! If FZEROI(NW) = 1, then the k-coefficients are zero - skip to the 110 112 ! special Gauss point at the end. 111 113 FZERO = FZEROI(NW) 112 114 113 115 IF(FZERO.ge.0.99) goto 40 114 116 115 117 DO NG=1,L_NGAUSS-1 116 118 117 119 if(TAUGSURF(NW,NG).lt. TLIMIT) then 118 120 fzero = fzero + (1.0D0-FZEROI(NW))*GWEIGHT(NG) 119 121 goto 30 120 122 end if 121 123 122 124 ! SET UP THE UPPER AND LOWER BOUNDARY CONDITIONS ON THE IR 123 125 ! CALCULATE THE DOWNWELLING RADIATION AT THE TOP OF THE MODEL 124 126 ! OR THE TOP LAYER WILL COOL TO SPACE UNPHYSICALLY 125 127 126 128 ! TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2)) 127 129 TAUTOP = TAUCUMI(2,NW,NG) 128 130 BTOP = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP 129 131 130 132 ! WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM 131 133 ! CALL A SUBROUTINE THAT SOLVES FOR THE FLUX TERMS 132 ! WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER 133 134 ! WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER 135 134 136 CALL GFLUXI(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG), 135 137 * TAUCUMI(1,NW,NG), 136 138 * WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP, 137 139 * BSURF,FTOPUP,FMUPI,FMDI) 138 140 139 141 ! NOW CALCULATE THE CUMULATIVE IR NET FLUX 140 142 NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*GWEIGHT(NG) 141 143 * * (1.0D0-FZEROI(NW)) 142 144 143 145 ! and same thing by spectral band... (RDW) 144 146 NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW) + FTOPUP * DWNI(NW) 145 147 * * GWEIGHT(NG) * (1.0D0-FZEROI(NW)) 146 148 147 149 DO L=1,L_NLEVRAD-1 148 150 ! CORRECT FOR THE WAVENUMBER INTERVALS 149 151 FMNETI(L) = FMNETI(L) + (FMUPI(L)-FMDI(L)) * DWNI(NW) 150 152 * * GWEIGHT(NG)*(1.0D0-FZEROI(NW)) 153 FMNETI_NU(L,NW) = FMNETI_NU(L,NW) + (FMUPI(L)-FMDI(L)) 154 * * DWNI(NW) * GWEIGHT(NG)*(1.0D0-FZEROI(NW)) 151 155 FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)*DWNI(NW)*GWEIGHT(NG) 152 156 * * (1.0D0-FZEROI(NW)) … … 157 161 * * GWEIGHT(NG) * (1.0D0 - FZEROI(NW)) 158 162 END DO 159 163 160 164 30 CONTINUE 161 165 162 166 END DO !End NGAUSS LOOP 163 167 164 168 40 CONTINUE 165 169 166 170 ! SPECIAL 17th Gauss point 167 171 NG = L_NGAUSS 168 172 169 173 ! TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2)) 170 174 TAUTOP = TAUCUMI(2,NW,NG) 171 175 BTOP = (1.0D0-EXP(-TAUTOP/UBARI))*PLTOP 172 176 173 177 ! WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM 174 178 ! CALL A SUBROUTINE THAT SOLVES FOR THE FLUX TERMS 175 ! WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER 176 179 ! WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER 180 177 181 CALL GFLUXI(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG), 178 182 * TAUCUMI(1,NW,NG), 179 183 * WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP, 180 184 * BSURF,FTOPUP,FMUPI,FMDI) 181 185 182 186 ! NOW CALCULATE THE CUMULATIVE IR NET FLUX 183 187 NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*FZERO 184 188 185 189 ! and same thing by spectral band... (RW) 186 190 NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW) 187 191 * +FTOPUP*DWNI(NW)*FZERO 188 192 189 193 DO L=1,L_NLEVRAD-1 190 194 ! CORRECT FOR THE WAVENUMBER INTERVALS 191 195 FMNETI(L) = FMNETI(L)+(FMUPI(L)-FMDI(L))*DWNI(NW)*FZERO 196 FMNETI_nu(L,NW) = FMNETI_nu(L,NW)+(FMUPI(L)-FMDI(L)) 197 * *DWNI(NW)*FZERO 192 198 FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)*DWNI(NW)*FZERO 193 199 FLUXDNI(L) = FLUXDNI(L) + FMDI(L)*DWNI(NW)*FZERO … … 196 202 * + FMUPI(L) * DWNI(NW) * FZERO 197 203 END DO 198 204 199 205 501 CONTINUE !End Spectral Interval LOOP 200 206 ! *** END OF MAJOR SPECTRAL INTERVAL LOOP IN THE INFRARED**** 201 207 202 208 END SUBROUTINE SFLUXI 203 209 -
trunk/LMDZ.PLUTO/libf/phypluto/sfluxv.F
r3184 r3275 2 2 3 3 implicit none 4 4 5 5 contains 6 6 7 7 SUBROUTINE SFLUXV(DTAUV,TAUV,TAUCUMV,RSFV,DWNV,WBARV,COSBV, 8 8 * UBAR0,STEL,NFLUXTOPV,FLUXTOPVDN, 9 * NFLUXOUTV_nu,NFLUXGNDV_nu, 10 * FMNETV,FLUXUPV,FLUXDNV,FZEROV,taugsurf) 9 * NFLUXOUTV_nu,NFLUXGNDV_nu,NFLUXTOPV_nu, 10 * FMNETV,FMNETV_NU,FLUXUPV,FLUXDNV, 11 * FZEROV,taugsurf) 11 12 12 13 use radinc_h, only: L_TAUMAX, L_LEVELS, L_NSPECTV, L_NGAUSS … … 18 19 19 20 real*8 FMNETV(L_NLAYRAD) 21 real*8 FMNETV_NU(L_NLAYRAD,L_NSPECTV) 20 22 real*8 TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS) 21 23 real*8 TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS) … … 27 29 real*8 FLUXUPV(L_NLAYRAD), FLUXDNV(L_NLAYRAD) 28 30 real*8 NFLUXTOPV, FLUXUP, FLUXDN,FLUXTOPVDN 31 real*8 NFLUXTOPV_nu(L_NSPECTV) 29 32 real*8 NFLUXOUTV_nu(L_NSPECTV) 30 33 real*8 NFLUXGNDV_nu(L_NSPECTV) … … 50 53 NFLUXOUTV_nu(NW)=0.0 51 54 NFLUXGNDV_nu(NW)=0.0 55 NFLUXTOPV_nu(nw)=0.0 56 DO L=1,L_NLAYRAD 57 FMNETV_NU(L,NW) = 0.0 58 END DO 52 59 END DO 53 60 … … 64 71 65 72 DO 500 NW=1,L_NSPECTV 66 73 67 74 F0PI = STEL(NW) 68 75 … … 84 91 BSURF = 0. 85 92 C LOOP OVER THE NTERMS BEGINNING HERE 86 93 87 94 88 95 ! FACTOR = 1.0D0 - WDEL(1)*CDEL(1)**2 … … 96 103 C CALL A SUBROUTINE THAT SOLVES FOR THE FLUX TERMS 97 104 C WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER 98 C 99 C FUW AND FDW ARE WORKING FLUX ARRAYS THAT WILL BE USED TO 105 C 106 C FUW AND FDW ARE WORKING FLUX ARRAYS THAT WILL BE USED TO 100 107 C RETURN FLUXES FOR A GIVEN NT 101 108 102 109 103 110 CALL GFLUXV(DTAUV(1,NW,NG),TAUV(1,NW,NG),TAUCUMV(1,NW,NG), 104 * WBARV(1,NW,NG),COSBV(1,NW,NG),UBAR0,F0PI,RSFV(NW), 111 * WBARV(1,NW,NG),COSBV(1,NW,NG),UBAR0,F0PI,RSFV(NW), 105 112 * BTOP,BSURF,FMUPV,FMDV,DIFFV,FLUXUP,FLUXDN) 106 113 107 C NOW CALCULATE THE CUMULATIVE VISIBLE NET FLUX 114 C NOW CALCULATE THE CUMULATIVE VISIBLE NET FLUX 108 115 109 116 NFLUXTOPV = NFLUXTOPV+(FLUXUP-FLUXDN)*GWEIGHT(NG)* … … 114 121 FMNETV(L)=FMNETV(L)+( FMUPV(L)-FMDV(L) )* 115 122 * GWEIGHT(NG)*(1.0-FZEROV(NW)) 123 FMNETV_NU(L,NW)=FMNETV_NU(L,NW)+( FMUPV(L)-FMDV(L) )* 124 * GWEIGHT(NG)*(1.0-FZEROV(NW)) 116 125 FLUXUPV(L) = FLUXUPV(L) + FMUPV(L)*GWEIGHT(NG)* 117 126 * (1.0-FZEROV(NW)) … … 120 129 END DO 121 130 131 c and same thing by spectral band... (RDW) 132 NFLUXTOPV_nu(NW) = NFLUXTOPV_nu(NW) 133 * +(FLUXUP-FLUXDN)*GWEIGHT(NG)* 134 * (1.0-FZEROV(NW)) 135 122 136 c band-resolved flux leaving TOA (RDW) 123 137 NFLUXOUTV_nu(NW) = NFLUXOUTV_nu(NW) … … 133 147 DIFFVT = DIFFVT + DIFFV*GWEIGHT(NG)*(1.0-FZEROV(NW)) 134 148 135 30 CONTINUE 136 137 END DO ! the Gauss loop 138 139 40 continue 149 30 CONTINUE 150 151 END DO ! the Gauss loop 152 153 40 continue 140 154 C Special 17th Gauss point 141 155 … … 143 157 144 158 C SET UP THE UPPER AND LOWER BOUNDARY CONDITIONS ON THE VISIBLE 145 159 146 160 BTOP = 0.0 147 161 148 162 C LOOP OVER THE NTERMS BEGINNING HERE 149 163 150 164 ETERM = MIN(TAUV(L_NLEVRAD,NW,NG),TAUMAX) 151 165 BSURF = RSFV(NW)*UBAR0*STEL(NW)*EXP(-ETERM/UBAR0) … … 155 169 C CALL A SUBROUTINE THAT SOLVES FOR THE FLUX TERMS 156 170 C WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER 157 C 158 C FUW AND FDW ARE WORKING FLUX ARRAYS THAT WILL BE USED TO 171 C 172 C FUW AND FDW ARE WORKING FLUX ARRAYS THAT WILL BE USED TO 159 173 C RETURN FLUXES FOR A GIVEN NT 160 174 … … 164 178 165 179 166 C NOW CALCULATE THE CUMULATIVE VISIBLE NET FLUX 180 C NOW CALCULATE THE CUMULATIVE VISIBLE NET FLUX 167 181 168 182 NFLUXTOPV = NFLUXTOPV+(FLUXUP-FLUXDN)*FZERO … … 170 184 DO L=1,L_NLAYRAD 171 185 FMNETV(L)=FMNETV(L)+( FMUPV(L)-FMDV(L) )*FZERO 186 FMNETV_NU(L,NW)=FMNETV_NU(L,NW)+( FMUPV(L)-FMDV(L) )*FZERO 172 187 FLUXUPV(L) = FLUXUPV(L) + FMUPV(L)*FZERO 173 188 FLUXDNV(L) = FLUXDNV(L) + FMDV(L)*FZERO 174 189 END DO 175 190 191 c and same thing by spectral band... (RDW) 192 NFLUXTOPV_nu(NW) = NFLUXTOPV_nu(NW) 193 * +(FLUXUP-FLUXDN)*FZERO 194 176 195 c band-resolved flux leaving TOA (RDW) 177 196 NFLUXOUTV_nu(NW) = NFLUXOUTV_nu(NW) … … 196 215 197 216 end module sfluxv_mod 198 217 -
trunk/LMDZ.PLUTO/libf/phypluto/surfdat_h.F90
r3195 r3275 17 17 !$OMP THREADPRIVATE(zmea,zstd,zsig,zgam,zthe) 18 18 real ttop 19 real,allocatable,dimension(:) :: kp ! TB ref pressure 19 real,allocatable,dimension(:) :: kp ! TB ref pressure 20 20 real p00 21 21 !$OMP THREADPRIVATE(ttop,kp,p00) 22 22 ! surface properties ! TB16 23 23 real alb_n2b,alb_n2a,alb_ch4,alb_co,alb_tho,emis_n2b,emis_n2a 24 !$OMP THREADPRIVATE(alb_n2b,alb_n2a,alb_ch4,alb_co,alb_tho,emis_n2b,emis_n2a) 24 !$OMP THREADPRIVATE(alb_n2b,alb_n2a,alb_ch4,alb_co,alb_tho,emis_n2b,emis_n2a) 25 25 real emis_ch4,emis_co,emis_tho,emis_tho_eq,alb_tho_eq,alb_ch4_eq 26 !$OMP THREADPRIVATE(emis_ch4,emis_co,emis_tho,emis_tho_eq,alb_tho_eq,alb_ch4_eq) 26 !$OMP THREADPRIVATE(emis_ch4,emis_co,emis_tho,emis_tho_eq,alb_tho_eq,alb_ch4_eq) 27 27 real ITN2,ITCH4,ITH2O,ITN2d,ITCH4d,ITH2Od,alb_ch4_s,albspe,emispe 28 !$OMP THREADPRIVATE(ITN2,ITCH4,ITH2O,ITN2d,ITCH4d,ITH2Od,alb_ch4_s,albspe,emispe) 28 !$OMP THREADPRIVATE(ITN2,ITCH4,ITH2O,ITN2d,ITCH4d,ITH2Od,alb_ch4_s,albspe,emispe) 29 29 real alb_tho_spe,emis_tho_spe 30 !$OMP THREADPRIVATE(alb_tho_spe,emis_tho_spe) 30 !$OMP THREADPRIVATE(alb_tho_spe,emis_tho_spe) 31 31 32 32 real,allocatable,dimension(:) :: dryness !"Dryness coefficient" for grnd water ice sublimation -
trunk/LMDZ.PLUTO/libf/phypluto/tracer_h.F90
r3237 r3275 51 51 ! Lists of constants for condensable tracers 52 52 real, save, allocatable :: constants_mass(:) ! molecular mass of the specie (g/mol) 53 real, save, allocatable :: constants_delta_ vapH(:) ! Enthalpy of vaporization (J/mol)53 real, save, allocatable :: constants_delta_gasH(:) ! Enthalpy of vaporization (J/mol) 54 54 real, save, allocatable :: constants_Tref(:) ! Ref temperature for Clausis-Clapeyron (K) 55 55 real, save, allocatable :: constants_Pref(:) ! Reference pressure for Clausius Clapeyron (Pa) … … 58 58 real, save, allocatable :: constants_metallicity_coeff(:) ! Coefficient to take into account the metallicity 59 59 real, save, allocatable :: constants_RCPV_generic(:) ! specific heat capacity of the tracer vapor at Tref 60 !$OMP THREADPRIVATE(constants_mass,constants_delta_ vapH,constants_Tref)60 !$OMP THREADPRIVATE(constants_mass,constants_delta_gasH,constants_Tref) 61 61 !$OMP THREADPRIVATE(constants_Pref,constants_epsi_generic) 62 62 !$OMP THREADPRIVATE(constants_RLVTT_generic,constants_metallicity_coeff,constants_RCPV_generic) -
trunk/LMDZ.PLUTO/libf/phypluto/vdifc_mod.F
r3195 r3275 435 435 do iq=1,nq 436 436 437 ! if (iq.ne.igcm_h2o_ vap) then437 ! if (iq.ne.igcm_h2o_gas) then 438 438 439 439 DO ig=1,ngrid
Note: See TracChangeset
for help on using the changeset viewer.