subroutine SISVAT_TS2 ! #ES. (ETSo_0,ETSo_1,ETSo_d) ! +------------------------------------------------------------------------+ ! | MAR SISVAT_TS2 Mon 16-08-2009 MAR | ! | SubRoutine SISVAT_TS2 computes the Soil/Snow temperature and fluxes | ! | using the same method as in LMDZ for consistency. | ! | The corresponding LMDZ routines are soil (soil.F90) and calcul_fluxs | ! | (calcul_fluxs_mod.F90). | ! +------------------------------------------------------------------------+ ! | | ! | | ! | PARAMETERS: klonv: Total Number of columns = | ! | ^^^^^^^^^^ = Total Number of grid boxes of surface type | ! | (land ice for now) | ! | | ! | INPUT: isnoSV = total Nb of Ice/Snow Layers | ! | ^^^^^ sol_SV : Downward Solar Radiation [W/m2] | ! | IRd_SV : Surface Downward Longwave Radiation [W/m2] | ! | VV__SV : SBL Top Wind Speed [m/s] | ! | TaT_SV : SBL Top Temperature [K] | ! | QaT_SV : SBL Top Specific Humidity [kg/kg] | ! | dzsnSV : Snow Layer Thickness [m] | ! | dt__SV : Time Step [s] | ! | | ! | SoSosv : Absorbed Solar Radiation by Surfac.(Normaliz)[-] | ! | Eso_sv : Soil+Snow Emissivity [-] | ! | ? rah_sv : Aerodynamic Resistance for Heat [s/m] | ! | | ! | dz1_SV : "inverse" layer thickness (centered) [1/m] | ! | dz2_SV : layer thickness (layer above (?)) [m] | ! | AcoHSV : coefficient for Enthalpy evolution, from atm. | ! | AcoHSV : coefficient for Enthalpy evolution, from atm. | ! | AcoQSV : coefficient for Humidity evolution, from atm. | ! | BcoQSV : coefficient for Humidity evolution, from atm. | ! | ps__SV : surface pressure [Pa] | ! | p1l_SV : 1st atmospheric layer pressure [Pa] | ! | cdH_SV : drag coeff Energy (?) | ! | rsolSV : Radiation balance surface [W/m2] | ! | lambSV : Coefficient for soil layer geometry [-] | ! | | ! | INPUT / TsisSV : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)| ! | OUTPUT: & Snow Temperatures (layers 1,2,...,nsno) [K] | ! | ^^^^^^ rsolSV : Radiation balance surface [W/m2] | ! | | ! | OUTPUT: IRs_SV : Soil IR Radiation [W/m2] | ! | ^^^^^^ HSs_sv : Sensible Heat Flux [W/m2] | ! | HLs_sv : Latent Heat Flux [W/m2] | ! | TsfnSV : new surface temperature [K] | ! | Evp_sv : Evaporation [kg/m2] | ! | dSdTSV : Sensible Heat Flux temp. derivative [W/m2/K] | ! | dLdTSV : Latent Heat Flux temp. derivative [W/m2/K] | ! | | ! | ? ETSo_0 : Snow/Soil Energy Power, before Forcing [W/m2] | ! | ? ETSo_1 : Snow/Soil Energy Power, after Forcing [W/m2] | ! | ? ETSo_d : Snow/Soil Energy Power Forcing [W/m2] | ! | | ! |________________________________________________________________________| USE yoethf_mod_h USE VAR_SV USE VARdSV USE VARySV USE VARtSV USE VARxSV USE VARphy USE indice_sol_mod USE yomcst_mod_h USE comsoil_mod_h IMPLICIT NONE ! +--Global Variables ! + ================ INCLUDE "FCTTRE.h" ! +--OUTPUT for Stand Alone NetCDF File ! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! #NC real*8 SOsoKL(klonv) ! Absorbed Solar Radiation ! #NC real*8 IRsoKL(klonv) ! Absorbed IR Radiation ! #NC real*8 HSsoKL(klonv) ! Absorbed Sensible Heat Flux ! #NC real*8 HLsoKL(klonv) ! Absorbed Latent Heat Flux ! #NC real*8 HLs_KL(klonv) ! Evaporation ! #NC real*8 HLv_KL(klonv) ! Transpiration ! #NC common/DumpNC/SOsoKL,IRsoKL ! #NC . ,HSsoKL,HLsoKL ! #NC . ,HLs_KL,HLv_KL ! +--Internal Variables ! + ================== integer :: ig,jk,isl real :: mu real :: Tsrf(klonv) ! surface temperature as extrapolated from soil real :: mug(klonv) !hj coef top layers real :: ztherm_i(klonv),zdz2(klonv,-nsol:nsno),z1s real :: pfluxgrd(klonv), pcapcal(klonv), cal(klonv) real :: beta(klonv), dif_grnd(klonv) real :: C_coef(klonv,-nsol:nsno),D_coef(klonv,-nsol:nsno) REAL, DIMENSION(klonv) :: zx_mh, zx_nh, zx_oh REAL, DIMENSION(klonv) :: zx_mq, zx_nq, zx_oq REAL, DIMENSION(klonv) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef REAL, DIMENSION(klonv) :: zx_sl, zx_k1 REAL, DIMENSION(klonv) :: d_ts REAL :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh REAL :: qsat_new, q1_new ! REAL, PARAMETER :: t_grnd = 271.35, t_coup = 273.15 ! REAL, PARAMETER :: max_eau_sol = 150.0 REAL, DIMENSION(klonv) :: IRs__D, dIRsdT REAL :: t_grnd ! not used parameter(t_grnd = 271.35) ! REAL :: t_coup ! distinguish evap/sublimation parameter(t_coup = 273.15) ! REAL :: max_eau_sol parameter(max_eau_sol = 150.0) ! write(*,*)'T check' ! ! DO ig = 1,knonv ! DO jk = 1,isnoSV(ig) !nsno ! IF (TsisSV(ig,jk) <= 1.) THEN !hj check ! TsisSV(ig,jk) = TsisSV(ig,isnoSV(ig)) ! ENDIF ! ! IF (TsisSV(ig,jk) <= 1.) THEN !hj check ! TsisSV(ig,jk) = 273.15 ! ENDIF ! END DO ! END DO !!======================================================================= !! I. First part: corresponds to soil.F90 in LMDZ !!======================================================================= DO ig = 1,knonv DO jk =1,isnoSV(ig) dz2_SV(ig,jk)=dzsnSV(ig,jk) !! use arithmetic center between layers to derive dz1 for snow layers for simplicity: dz1_SV(ig,jk)=2./(dzsnSV(ig,jk)+dzsnSV(ig,jk-1)) ENDDO ENDDO DO ig = 1,knonv ztherm_i(ig) = inertie_lic IF (isnoSV(ig) > 0) ztherm_i(ig) = inertie_sno ENDDO !!----------------------------------------------------------------------- !! 1) !! Calculation of Cgrf and Dgrd coefficients using soil temperature from !! previous time step. !! !! These variables are recalculated on the local compressed grid instead !! of saved in restart file. !!----------------------------------------------------------------------- DO ig=1,knonv DO jk=-nsol,nsno zdz2(ig,jk)=dz2_SV(ig,jk)/dt__SV !ptimestep ENDDO ENDDO DO ig=1,knonv z1s = zdz2(ig,-nsol)+dz1_SV(ig,-nsol+1) C_coef(ig,-nsol+1)=zdz2(ig,-nsol)*TsisSV(ig,-nsol)/z1s D_coef(ig,-nsol+1)=dz1_SV(ig,-nsol+1)/z1s ENDDO DO ig=1,knonv DO jk=-nsol+1,isnoSV(ig)-1,1 z1s = 1./(zdz2(ig,jk)+dz1_SV(ig,jk+1)+dz1_SV(ig,jk) & *(1.-D_coef(ig,jk))) C_coef(ig,jk+1)= & (TsisSV(ig,jk)*zdz2(ig,jk) & +dz1_SV(ig,jk)*C_coef(ig,jk)) * z1s D_coef(ig,jk+1)=dz1_SV(ig,jk+1)*z1s ENDDO ENDDO !!----------------------------------------------------------------------- !! 2) !! Computation of the soil temperatures using the Cgrd and Dgrd !! coefficient computed above !! !!----------------------------------------------------------------------- !! Extrapolate surface Temperature !hj check mu=1./((2.**1.5-1.)/(2.**(0.5)-1.)-1.) ! IF (knonv>0) THEN ! DO ig=1,8 ! write(*,*)ig,'sisvat: Tsis ',TsisSV(ig,isnoSV(ig)) ! write(*,*)'max-1 ',TsisSV(ig,isnoSV(ig)-1) ! write(*,*)'max-2 ',TsisSV(ig,isnoSV(ig)-2) ! write(*,*)'0 ',TsisSV(ig,0) !! write(*,*)min(max(isnoSV(ig),0),1),max(1-isnoSV(ig),0) ! ENDDO ! END IF DO ig=1,knonv IF (isnoSV(ig).GT.0) THEN IF (isnoSV(ig).GT.1) THEN mug(ig)=1./(1.+dzsnSV(ig,isnoSV(ig)-1)/dzsnSV(ig,isnoSV(ig))) !mu ELSE mug(ig) = 1./(1.+dzsnSV(ig,isnoSV(ig)-1)/dz_dSV(0)) !mu ENDIF ELSE mug(ig) = lambSV ENDIF IF (mug(ig) .LE. 0.05) THEN write(*,*)'Attention mu low', mug(ig) ENDIF IF (mug(ig) .GE. 0.98) THEN write(*,*)'Attention mu high', mug(ig) ENDIF Tsrf(ig)=(1.5*TsisSV(ig,isnoSV(ig))-0.5*TsisSV(ig,isnoSV(ig)-1))& *min(max(isnoSV(ig),0),1)+ & ((mug(ig)+1)*TsisSV(ig,0)-mug(ig)*TsisSV(ig,-1)) & *max(1-isnoSV(ig),0) ENDDO !! Surface temperature DO ig=1,knonv TsisSV(ig,isnoSV(ig))=(mug(ig)*C_coef(ig,isnoSV(ig))+Tsf_SV(ig))/ & (mug(ig)*(1.-D_coef(ig,isnoSV(ig)))+1.) ENDDO !! Other temperatures DO ig=1,knonv DO jk=isnoSV(ig),-nsol+1,-1 TsisSV(ig,jk-1)=C_coef(ig,jk)+D_coef(ig,jk) & *TsisSV(ig,jk) ENDDO ENDDO ! write(*,*)ig,'Tsis',TsisSV(ig,0) ! IF (indice == is_sic) THEN ! DO ig = 1,knonv ! TsisSV(ig,-nsol) = RTT - 1.8 ! END DO ! ENDIF !C !hj new 11 03 2010 DO ig=1,knonv isl = isnoSV(ig) ! dIRsdT(ig) = Eso_sv(ig)* SteBo * 4. & ! - d(IR)/d(T) ! & * Tsf_SV(ig) & !T TsisSV(ig,isl) ! ! & * Tsf_SV(ig) & !TsisSV(ig,isl) ! ! & * Tsf_SV(ig) !TsisSV(ig,isl) ! ! IRs__D(ig) = dIRsdT(ig)* Tsf_SV(ig) * 0.75 !TsisSV(ig,isl) * 0.75 !: dIRsdT(ig) = Eso_sv(ig)* StefBo * 4. & ! - d(IR)/d(T) * TsisSV(ig,isl) & * TsisSV(ig,isl) & * TsisSV(ig,isl) IRs__D(ig) = dIRsdT(ig)* TsisSV(ig,isl) * 0.75 !: END DO ! !hj !!----------------------------------------------------------------------- !! 3) !! Calculate the Cgrd and Dgrd coefficient corresponding to actual soil !! temperature !!----------------------------------------------------------------------- DO ig=1,knonv z1s = zdz2(ig,-nsol)+dz1_SV(ig,-nsol+1) C_coef(ig,-nsol+1) = zdz2(ig,-nsol)*TsisSV(ig,-nsol)/z1s D_coef(ig,-nsol+1) = dz1_SV(ig,-nsol+1)/z1s ENDDO DO ig=1,knonv DO jk=-nsol+1,isnoSV(ig)-1,1 z1s = 1./(zdz2(ig,jk)+dz1_SV(ig,jk+1)+dz1_SV(ig,jk) & *(1.-D_coef(ig,jk))) C_coef(ig,jk+1) = (TsisSV(ig,jk)*zdz2(ig,jk)+ & dz1_SV(ig,jk)*C_coef(ig,jk)) * z1s D_coef(ig,jk+1) = dz1_SV(ig,jk+1)*z1s ENDDO ENDDO !!----------------------------------------------------------------------- !! 4) !! Computation of the surface diffusive flux from ground and !! calorific capacity of the ground !!----------------------------------------------------------------------- DO ig=1,knonv !! (pfluxgrd) pfluxgrd(ig) = ztherm_i(ig)*dz1_SV(ig,isnoSV(ig))* & (C_coef(ig,isnoSV(ig))+(D_coef(ig,isnoSV(ig))-1.) & *TsisSV(ig,isnoSV(ig))) !! (pcapcal) pcapcal(ig) = ztherm_i(ig)* & (dz2_SV(ig,isnoSV(ig))+dt__SV*(1.-D_coef(ig,isnoSV(ig))) & *dz1_SV(ig,isnoSV(ig))) z1s = mug(ig)*(1.-D_coef(ig,isnoSV(ig)))+1. pcapcal(ig) = pcapcal(ig)/z1s pfluxgrd(ig) = ( pfluxgrd(ig) & + pcapcal(ig) * (TsisSV(ig,isnoSV(ig)) * z1s & - mug(ig)* C_coef(ig,isnoSV(ig)) & - Tsf_SV(ig)) /dt__SV ) ENDDO cal(1:knonv) = RCPD / pcapcal(1:knonv) rsolSV(1:knonv) = rsolSV(1:knonv) + pfluxgrd(1:knonv) !!======================================================================= !! II. Second part: corresponds to calcul_fluxs_mod.F90 in LMDZ !!======================================================================= Evp_sv = 0. ! #NC HSsoKL=0. ! #NC HLsoKL=0. dSdTSV = 0. dLdTSV = 0. beta(:) = 1.0 dif_grnd(:) = 0.0 !! zx_qs = qsat en kg/kg !!**********************************************************************x*************** DO ig = 1,knonv IF (ps__SV(ig).LT.1.) THEN ! write(*,*)'ig',ig,'ps',ps__SV(ig) ps__SV(ig)=max(ps__SV(ig),1.e-8) ENDIF IF (p1l_SV(ig).LT.1.) THEN ! write(*,*)'ig',ig,'p1l',p1l_SV(ig) p1l_SV(ig)=max(p1l_SV(ig),1.e-8) ENDIF IF (TaT_SV(ig).LT.180.) THEN ! write(*,*)'ig',ig,'TaT',TaT_SV(ig) TaT_SV(ig)=max(TaT_SV(ig),180.) ENDIF IF (QaT_SV(ig).LT.1.e-8) THEN ! write(*,*)'ig',ig,'QaT',QaT_SV(ig) QaT_SV(ig)=max(QaT_SV(ig),1.e-8) ENDIF IF (Tsf_SV(ig).LT.100.) THEN ! write(*,*)'ig',ig,'Tsf',Tsf_SV(ig) Tsf_SV(ig)=max(Tsf_SV(ig),180.) ENDIF IF (Tsf_SV(ig).GT.500.) THEN ! write(*,*)'ig',ig,'Tsf',Tsf_SV(ig) Tsf_SV(ig)=min(Tsf_SV(ig),400.) ENDIF ! IF (Tsrf(ig).LT.1.) THEN !! write(*,*)'ig',ig,'Tsrf',Tsrf(ig) ! Tsrf(ig)=max(Tsrf(ig),TaT_SV(ig)-20.) ! ENDIF IF (cdH_SV(ig).LT.1.e-10) THEN ! IF (ig.le.3) write(*,*)'ig',ig,'cdH',cdH_SV(ig) cdH_SV(ig)=.5 ENDIF ENDDO DO ig = 1,knonv zx_pkh(ig) = 1. ! (ps__SV(ig)/ps__SV(ig))**RKAPPA IF (thermcep) THEN zdelta=MAX(0.,SIGN(1.,rtt-Tsf_SV(ig))) zcvm5 = R5LES*LhvH2O*(1.-zdelta) + R5IES*LhsH2O*zdelta zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*QaT_SV(ig)) zx_qs= r2es * FOEEW(Tsf_SV(ig),zdelta)/ps__SV(ig) zx_qs=MIN(0.5,zx_qs) ! !write(*,*)'zcor',retv*zx_qs zcor=1./(1.-retv*zx_qs) zx_qs=zx_qs*zcor zx_dq_s_dh = FOEDE(Tsf_SV(ig),zdelta,zcvm5,zx_qs,zcor) & /LhvH2O / zx_pkh(ig) ELSE IF (Tsf_SV(ig).LT.t_coup) THEN zx_qs = qsats(Tsf_SV(ig)) / ps__SV(ig) zx_dq_s_dh = dqsats(Tsf_SV(ig),zx_qs)/LhvH2O & / zx_pkh(ig) ELSE zx_qs = qsatl(Tsf_SV(ig)) / ps__SV(ig) zx_dq_s_dh = dqsatl(Tsf_SV(ig),zx_qs)/LhvH2O & / zx_pkh(ig) ENDIF ENDIF zx_dq_s_dt(ig) = RCPD * zx_pkh(ig) * zx_dq_s_dh zx_qsat(ig) = zx_qs ! zx_coef(ig) = cdH_SV(ig) * & ! & (1.0+SQRT(u1lay(ig)**2+v1lay(ig)**2)) * & ! & p1l_SV(ig)/(RD*t1lay(ig)) zx_coef(ig) = cdH_SV(ig) * & (1.0+VV__SV(ig)) * & p1l_SV(ig)/(RD*TaT_SV(ig)) ENDDO !! === Calcul de la temperature de surface === !! zx_sl = chaleur latente d'evaporation ou de sublimation !!**************************************************************************** DO ig = 1,knonv zx_sl(ig) = LhvH2O IF (Tsf_SV(ig) .LT. RTT) zx_sl(ig) = LhsH2O zx_k1(ig) = zx_coef(ig) ENDDO DO ig = 1,knonv !! Q zx_oq(ig) = 1. - (beta(ig) * zx_k1(ig) * BcoQSV(ig) * dt__SV) zx_mq(ig) = beta(ig) * zx_k1(ig) * & (AcoQSV(ig) - zx_qsat(ig) + & zx_dq_s_dt(ig) * Tsf_SV(ig)) & / zx_oq(ig) zx_nq(ig) = beta(ig) * zx_k1(ig) * (-1. * zx_dq_s_dt(ig)) & / zx_oq(ig) !! H zx_oh(ig) = 1. - (zx_k1(ig) * BcoHSV(ig) * dt__SV) zx_mh(ig) = zx_k1(ig) * AcoHSV(ig) / zx_oh(ig) zx_nh(ig) = - (zx_k1(ig) * RCPD * zx_pkh(ig))/ zx_oh(ig) !! surface temperature TsfnSV(ig) = (Tsf_SV(ig) + cal(ig)/RCPD * zx_pkh(ig) * dt__SV * & (rsolSV(ig) + zx_mh(ig) + zx_sl(ig) * zx_mq(ig)) & + dif_grnd(ig) * t_grnd * dt__SV)/ & ( 1. - dt__SV * cal(ig)/(RCPD * zx_pkh(ig)) * & (zx_nh(ig) + zx_sl(ig) * zx_nq(ig)) & + dt__SV * dif_grnd(ig)) !hj rajoute 22 11 2010 tuning... TsfnSV(ig) = min(RTT+0.02,TsfnSV(ig)) d_ts(ig) = TsfnSV(ig) - Tsf_SV(ig) !!== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas !!== flux_t est le flux de cpt (energie sensible): j/(m**2 s) Evp_sv(ig) = - zx_mq(ig) - zx_nq(ig) * TsfnSV(ig) HLs_sv(ig) = - Evp_sv(ig) * zx_sl(ig) HSs_sv(ig) = zx_mh(ig) + zx_nh(ig) * TsfnSV(ig) !! Derives des flux dF/dTs (W m-2 K-1): dSdTSV(ig) = zx_nh(ig) dLdTSV(ig) = zx_sl(ig) * zx_nq(ig) !hj new 11 03 2010 isl = isnoSV(ig) ! TsisSV(ig,isl) = TsfnSV(ig) IRs_SV(ig) = IRs__D(ig) & ! - dIRsdT(ig) * TsfnSV(ig) !TsisSV(ig,isl)? ! ! hj ! #NC SOsoKL(ig) = sol_SV(ig) * SoSosv(ig) ! Absorbed Sol. ! #NC IRsoKL(ig) = IRs_SV(ig) & !Up Surf. IR ! #NC& + tau_sv(ig) *IRd_SV(ig)*Eso_sv(ig) & !Down Atm IR ! #NC& -(1.0-tau_sv(ig)) *0.5*IRv_sv(ig) ! Down Veg IR ! #NC HLsoKL(ig) = HLs_sv(ig) ! #NC HSsoKL(ig) = HSs_sv(ig) ! #NC HLs_KL(ig) = Evp_sv(ig) !! Nouvelle valeure de l'humidite au dessus du sol qsat_new=zx_qsat(ig) + zx_dq_s_dt(ig) * d_ts(ig) q1_new = AcoQSV(ig) - BcoQSV(ig)* Evp_sv(ig)*dt__SV QaT_SV(ig)=q1_new*(1.-beta(ig)) + beta(ig)*qsat_new ENDDO end subroutine sisvat_ts2