Index: LMDZ6/trunk/libf/phylmdiso/cdrag.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/cdrag.F90 (revision 3940) +++ (revision ) @@ -1,411 +1,0 @@ -! -!$Id: cdrag.F90 2868 2017-05-02 11:16:06Z fhourdin $ -! - SUBROUTINE cdrag(knon, nsrf, & - speed, t1, q1, zgeop1, & - psol, tsurf, qsurf, z0m, z0h, & - cdm, cdh, zri, pref) - - USE dimphy - USE indice_sol_mod - USE print_control_mod, ONLY: lunout, prt_level - USE ioipsl_getin_p_mod, ONLY : getin_p - - IMPLICIT NONE -! ================================================================= c -! -! Objet : calcul des cdrags pour le moment (pcfm) et -! les flux de chaleur sensible et latente (pcfh) d'apr??s -! Louis 1982, Louis 1979, King et al 2001 -! ou Zilitinkevich et al 2002 pour les cas stables, Louis 1979 -! et 1982 pour les cas instables -! -! Modified history: -! writting on the 20/05/2016 -! modified on the 13/12/2016 to be adapted to LMDZ6 -! -! References: -! Louis, J. F., 1979: A parametric model of vertical eddy fluxes in the -! atmosphere. Boundary-Layer Meteorology. 01/1979; 17(2):187-202. -! Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982: `A short history of the -! operational PBL parametrization at ECMWF'. Workshop on boundary layer -! parametrization, November 1981, ECMWF, Reading, England. -! Page: 19. Equations in Table 1. -! Mascart P, Noilhan J, Giordani H 1995.A MODIFIED PARAMETERIZATION OF FLUX-PROFILE RELATIONSHIPS -! IN THE SURFACE LAYER USING DIFFERENT ROUGHNESS LENGTH VALUES FOR HEAT AND MOMENTUM -! Boundary-Layer Meteorology 72: 331-344 -! Anton Beljaars. May 1992. The parametrization of the planetary boundary layer. -! European Centre for Medium-Range Weather Forecasts. -! Equations: 110-113. Page 40. -! Miller,M.J., A.C.M.Beljaars, T.N.Palmer. 1992. The sensitivity of the ECMWF -! model to the parameterization of evaporation from the tropical oceans. J. -! Climate, 5:418-434. -! King J.C, Connolley, W.M ad Derbyshire S.H. 2001, Sensitivity of Modelled Antarctic climate -! to surface and boundary-layer flux parametrizations -! QJRMS, 127, pp 779-794 -! -! ================================================================= c -! ================================================================= c -! On choisit le couple de fonctions de correction avec deux flags: -! Un pour les cas instables, un autre pour les cas stables -! -! iflag_corr_insta: -! 1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h) -! 2: Louis 1982 -! 3: Laurent Li -! -! iflag_corr_sta: -! 1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h) -! 2: Louis 1982 -! 3: Laurent Li -! 4: King 2001 (SHARP) -! 5: MO 1st order theory (allow collapse of turbulence) -! -! -!***************************************************************** -! Parametres d'entree -!***************************************************************** - - INTEGER, INTENT(IN) :: knon, nsrf ! nombre de points de grille sur l'horizontal + type de surface - REAL, DIMENSION(klon), INTENT(IN) :: speed ! module du vent au 1er niveau du modele - REAL, DIMENSION(klon), INTENT(IN) :: zgeop1! geopotentiel au 1er niveau du modele - REAL, DIMENSION(klon), INTENT(IN) :: tsurf ! Surface temperature (K) - REAL, DIMENSION(klon), INTENT(IN) :: qsurf ! Surface humidity (Kg/Kg) - REAL, DIMENSION(klon), INTENT(IN) :: z0m, z0h ! Rugosity at surface (m) - REAL, DIMENSION(klon), INTENT(IN) :: t1 ! Temperature au premier niveau (K) - REAL, DIMENSION(klon), INTENT(IN) :: q1 ! humidite specifique au premier niveau (kg/kg) - REAL, DIMENSION(klon), INTENT(IN) :: psol ! pression au sol - - - -! Parametres de sortie -!****************************************************************** - REAL, DIMENSION(klon), INTENT(OUT) :: cdm ! Drag coefficient for heat flux - REAL, DIMENSION(klon), INTENT(OUT) :: cdh ! Drag coefficient for momentum - REAL, DIMENSION(klon), INTENT(OUT) :: zri ! Richardson number - REAL, DIMENSION(klon), INTENT(OUT) :: pref ! Pression au niveau zgeop/RG - -! Variables Locales -!****************************************************************** - - - INCLUDE "YOMCST.h" - INCLUDE "YOETHF.h" - INCLUDE "clesphys.h" - - - REAL, PARAMETER :: CKAP=0.40, CKAPT=0.42 - REAL CEPDU2 - REAL ALPHA - REAL CB,CC,CD,C2,C3 - REAL MU, CM, CH, B, CMstar, CHstar - REAL PM, PH, BPRIME - REAL C - INTEGER ng_q1 ! Number of grids that q1 < 0.0 - INTEGER ng_qsurf ! Number of grids that qsurf < 0.0 - INTEGER i - REAL zdu2, ztsolv - REAL ztvd, zscf - REAL zucf, zcr - REAL friv, frih - REAL, DIMENSION(klon) :: FM, FH ! stability functions - REAL, DIMENSION(klon) :: cdmn, cdhn ! Drag coefficient in neutral conditions - REAL zzzcd - - LOGICAL, SAVE :: firstcall = .TRUE. - !$OMP THREADPRIVATE(firstcall) - INTEGER, SAVE :: iflag_corr_sta - !$OMP THREADPRIVATE(iflag_corr_sta) - INTEGER, SAVE :: iflag_corr_insta - !$OMP THREADPRIVATE(iflag_corr_insta) - -!===================================================================c -! Valeurs numeriques des constantes -!===================================================================c - - -! Minimum du carre du vent - - CEPDU2 = (0.1)**2 - -! Louis 1982 - - CB=5.0 - CC=5.0 - CD=5.0 - - -! King 2001 - - C2=0.25 - C3=0.0625 - - -! Louis 1979 - - BPRIME=4.7 - B=9.4 - - -!MO - - ALPHA=5.0 - - -! ================================================================= c -! Tests avant de commencer -! Fuxing WANG, 04/03/2015 -! To check if there are negative q1, qsurf values. -!====================================================================c - ng_q1 = 0 ! Initialization - ng_qsurf = 0 ! Initialization - DO i = 1, knon - IF (q1(i).LT.0.0) ng_q1 = ng_q1 + 1 - IF (qsurf(i).LT.0.0) ng_qsurf = ng_qsurf + 1 - ENDDO - IF (ng_q1.GT.0 .and. prt_level > 5) THEN - WRITE(lunout,*)" *** Warning: Negative q1(humidity at 1st level) values in cdrag.F90 !" - WRITE(lunout,*)" The total number of the grids is: ", ng_q1 - WRITE(lunout,*)" The negative q1 is set to zero " -! abort_message="voir ci-dessus" -! CALL abort_physic(modname,abort_message,1) - ENDIF - IF (ng_qsurf.GT.0 .and. prt_level > 5) THEN - WRITE(lunout,*)" *** Warning: Negative qsurf(humidity at surface) values in cdrag.F90 !" - WRITE(lunout,*)" The total number of the grids is: ", ng_qsurf - WRITE(lunout,*)" The negative qsurf is set to zero " -! abort_message="voir ci-dessus" -! CALL abort_physic(modname,abort_message,1) - ENDIF - - - -!=============================================================================c -! Calcul du cdrag -!=============================================================================c - -! On choisit les fonctions de stabilite utilisees au premier appel -!************************************************************************** - IF (firstcall) THEN - iflag_corr_sta=2 - iflag_corr_insta=2 - - CALL getin_p('iflag_corr_sta',iflag_corr_sta) - CALL getin_p('iflag_corr_insta',iflag_corr_insta) - - firstcall = .FALSE. - ENDIF - -!xxxxxxxxxxxxxxxxxxxxxxx - DO i = 1, knon ! Boucle sur l'horizontal -!xxxxxxxxxxxxxxxxxxxxxxx - - -! calculs preliminaires: -!*********************** - - - zdu2 = MAX(CEPDU2, speed(i)**2) - pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* & - (1.+ RETV * max(q1(i),0.0)))) ! negative q1 set to zero - ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0)) ! negative qsurf set to zero - ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) & - *(1.+RETV*max(q1(i),0.0)) ! negative q1 set to zero - zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd) - - -! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)): -!******************************************************************** - - zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i))) - cdmn(i) = zzzcd*zzzcd - cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i)))) - - -! Calcul des fonctions de stabilit?? FMs, FHs, FMi, FHi : -!******************************************************* - -!'''''''''''''' -! Cas instables -!'''''''''''''' - - IF (zri(i) .LT. 0.) THEN - - - SELECT CASE (iflag_corr_insta) - - CASE (1) ! Louis 1979 + Mascart 1995 - - MU=LOG(MAX(z0m(i)/z0h(i),0.01)) - CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3) - PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3) - CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3) - PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3) - CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & - & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i))) & - & * ((zgeop1(i)/(RG*z0h(i)))**PH) - CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & - & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & - & * ((zgeop1(i)/(RG*z0m(i)))**PM) - - - - - FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i)))) - FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i)))) - - CASE (2) ! Louis 1982 - - zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) & - *(1.0+zgeop1(i)/(RG*z0m(i))))) - FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min) - FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min) - - - CASE (3) ! Laurent Li - - - FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min) - FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min) - - - - CASE default ! Louis 1982 - - zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) & - *(1.0+zgeop1(i)/(RG*z0m(i))))) - FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min) - FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min) - - - END SELECT - - - -! Calcul des drags - - - cdm(i)=cdmn(i)*FM(i) - cdh(i)=f_cdrag_ter*cdhn(i)*FH(i) - - -! Traitement particulier des cas oceaniques -! on applique Miller et al 1992 en l'absence de gustiness - - IF (nsrf == is_oce) THEN -! cdh(i)=f_cdrag_oce*cdhn(i)*FH(i) - - IF(iflag_gusts==0) THEN - zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.) - cdh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25) - ENDIF - - - cdm(i)=MIN(cdm(i),cdmmax) - cdh(i)=MIN(cdh(i),cdhmax) - - END IF - - - - ELSE - -!''''''''''''''' -! Cas stables : -!''''''''''''''' - zri(i) = MIN(20.,zri(i)) - - SELECT CASE (iflag_corr_sta) - - CASE (1) ! Louis 1979 + Mascart 1995 - - FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min) - FH(i)=FM(i) - - - CASE (2) ! Louis 1982 - - zscf = SQRT(1.+CD*ABS(zri(i))) - FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min) - FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min ) - - - CASE (3) ! Laurent Li - - FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min) - FH(i)=FM(i) - - - CASE (4) ! King 2001 - - if (zri(i) .LT. C2/2.) then - FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min) - FH(i)= FM(i) - - - else - FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min) - FH(i)= FM(i) - endif - - - CASE (5) ! MO - - if (zri(i) .LT. 1./alpha) then - - FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min) - FH(i)=FM(i) - - else - - - FM(i)=MAX(1E-7,f_ri_cd_min) - FH(i)=FM(i) - - endif - - - - - - CASE default ! Louis 1982 - - zscf = SQRT(1.+CD*ABS(zri(i))) - FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min) - FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min ) - - - - END SELECT - -! Calcul des drags - - - cdm(i)=cdmn(i)*FM(i) - cdh(i)=f_cdrag_ter*cdhn(i)*FH(i) - - IF(nsrf.EQ.is_oce) THEN - - cdh(i)=f_cdrag_oce*cdhn(i)*FH(i) - cdm(i)=MIN(cdm(i),cdmmax) - cdh(i)=MIN(cdh(i),cdhmax) - - ENDIF - - - ENDIF - - - - -!xxxxxxxxxxx - END DO ! Fin de la boucle sur l'horizontal -!xxxxxxxxxxx -! ================================================================= c - - - -END SUBROUTINE cdrag - - -!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Index: LMDZ6/trunk/libf/phylmdiso/hbtm.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/hbtm.F90 (revision 3940) +++ (revision ) @@ -1,758 +1,0 @@ - -! $Header$ - - -SUBROUTINE hbtm(knon, paprs, pplay, t2m, t10m, q2m, q10m, ustar, wstar, & - flux_t, flux_q, u, v, t, q, pblh, cape, eauliq, ctei, pblt, therm, trmb1, & - trmb2, trmb3, plcl) - USE dimphy - IMPLICIT NONE - - ! *************************************************************** - ! * * - ! * HBTM2 D'apres Holstag&Boville et Troen&Mahrt * - ! * JAS 47 BLM * - ! * Algorithme These Anne Mathieu * - ! * Critere d'Entrainement Peter Duynkerke (JAS 50) * - ! * written by : Anne MATHIEU & Alain LAHELLEC, 22/11/99 * - ! * features : implem. exces Mathieu * - ! *************************************************************** - ! * mods : decembre 99 passage th a niveau plus bas. voir fixer * - ! * la prise du th a z/Lambda = -.2 (max Ray) * - ! * Autre algo : entrainement ~ Theta+v =cste mais comment=>The?* - ! * on peut fixer q a .7qsat(cf non adiab)=>T2 et The2 * - ! * voir aussi //KE pblh = niveau The_e ou l = env. * - ! *************************************************************** - ! * fin therm a la HBTM passage a forme Mathieu 12/09/2001 * - ! *************************************************************** - ! * - - - ! AM Fev 2003 - ! Adaptation a LMDZ version couplee - - ! Pour le moment on fait passer en argument les grdeurs de surface : - ! flux, t,q2m, t,q10m, on va utiliser systematiquement les grdeurs a 2m ms - ! on garde la possibilite de changer si besoin est (jusqu'a present la - ! forme de HB avec le 1er niveau modele etait conservee) - - - - - - include "YOMCST.h" - REAL rlvcp, reps - ! Arguments: - - INTEGER knon ! nombre de points a calculer - ! AM - REAL t2m(klon), t10m(klon) ! temperature a 2 et 10m - REAL q2m(klon), q10m(klon) ! q a 2 et 10m - REAL ustar(klon) - REAL wstar(klon) ! w*, convective velocity scale - REAL paprs(klon, klev+1) ! pression a inter-couche (Pa) - REAL pplay(klon, klev) ! pression au milieu de couche (Pa) - REAL flux_t(klon, klev), flux_q(klon, klev) ! Flux - REAL u(klon, klev) ! vitesse U (m/s) - REAL v(klon, klev) ! vitesse V (m/s) - REAL t(klon, klev) ! temperature (K) - REAL q(klon, klev) ! vapeur d'eau (kg/kg) - ! AM REAL cd_h(klon) ! coefficient de friction au sol pour chaleur - ! AM REAL cd_m(klon) ! coefficient de friction au sol pour vitesse - - INTEGER isommet - ! um PARAMETER (isommet=klev) ! limite max sommet pbl - REAL, PARAMETER :: vk = 0.35 ! Von Karman => passer a .41 ! cf U.Olgstrom - REAL, PARAMETER :: ricr = 0.4 - REAL, PARAMETER :: fak = 8.5 ! b calcul du Prandtl et de dTetas - REAL, PARAMETER :: fakn = 7.2 ! a - REAL, PARAMETER :: onet = 1.0/3.0 - REAL, PARAMETER :: t_coup = 273.15 - REAL, PARAMETER :: zkmin = 0.01 - REAL, PARAMETER :: betam = 15.0 ! pour Phim / h dans la S.L stable - REAL, PARAMETER :: betah = 15.0 - REAL, PARAMETER :: betas = 5.0 ! Phit dans la S.L. stable (mais 2 formes / z/OBL<>1 - REAL, PARAMETER :: sffrac = 0.1 ! S.L. = z/h < .1 - REAL, PARAMETER :: usmin = 1.E-12 - REAL, PARAMETER :: binm = betam*sffrac - REAL, PARAMETER :: binh = betah*sffrac - REAL, PARAMETER :: ccon = fak*sffrac*vk - REAL, PARAMETER :: b1 = 70., b2 = 20. - REAL, PARAMETER :: zref = 2. ! Niveau de ref a 2m peut eventuellement - ! etre choisi a 10m - REAL q_star, t_star - REAL b212, b2sr ! Lambert correlations T' q' avec T* q* - - REAL z(klon, klev) - ! AM REAL pcfm(klon,klev), pcfh(klon,klev) - INTEGER i, k, j - REAL zxt - ! AM REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy - ! AM REAL zx_alf1, zx_alf2 ! parametres pour extrapolation - REAL khfs(klon) ! surface kinematic heat flux [mK/s] - REAL kqfs(klon) ! sfc kinematic constituent flux [m/s] - REAL heatv(klon) ! surface virtual heat flux - REAL rhino(klon, klev) ! bulk Richardon no. mais en Theta_v - LOGICAL unstbl(klon) ! pts w/unstbl pbl (positive virtual ht flx) - LOGICAL stblev(klon) ! stable pbl with levels within pbl - LOGICAL unslev(klon) ! unstbl pbl with levels within pbl - LOGICAL unssrf(klon) ! unstb pbl w/lvls within srf pbl lyr - LOGICAL unsout(klon) ! unstb pbl w/lvls in outer pbl lyr - LOGICAL check(klon) ! True=>chk if Richardson no.>critcal - LOGICAL omegafl(klon) ! flag de prolongerment cape pour pt Omega - REAL pblh(klon) - REAL pblt(klon) - REAL plcl(klon) - ! AM REAL cgh(klon,2:klev) ! counter-gradient term for heat [K/m] - ! AM REAL cgq(klon,2:klev) ! counter-gradient term for constituents - ! AM REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux) - REAL unsobklen(klon) ! Monin-Obukhov lengh - ! AM REAL ztvd, ztvu, - REAL zdu2 - REAL therm(klon) ! thermal virtual temperature excess - REAL trmb1(klon), trmb2(klon), trmb3(klon) - ! Algorithme thermique - REAL s(klon, klev) ! [P/Po]^Kappa milieux couches - REAL th_th(klon) ! potential temperature of thermal - REAL the_th(klon) ! equivalent potential temperature of thermal - REAL qt_th(klon) ! total water of thermal - REAL tbef(klon) ! T thermique niveau precedent - REAL qsatbef(klon) - LOGICAL zsat(klon) ! le thermique est sature - REAL cape(klon) ! Cape du thermique - REAL kape(klon) ! Cape locale - REAL eauliq(klon) ! Eau liqu integr du thermique - REAL ctei(klon) ! Critere d'instab d'entrainmt des nuages de CL - REAL the1, the2, aa, bb, zthvd, zthvu, xintpos, qqsat - ! IM 091204 BEG - REAL a1, a2, a3 - ! IM 091204 END - REAL xhis, rnum, denom, th1, th2, thv1, thv2, ql2 - REAL dqsat_dt, qsat2, qt1, q2, t1, t2, xnull, delt_the - REAL delt_qt, delt_2, quadsat, spblh, reduc - - REAL phiminv(klon) ! inverse phi function for momentum - REAL phihinv(klon) ! inverse phi function for heat - REAL wm(klon) ! turbulent velocity scale for momentum - REAL fak1(klon) ! k*ustar*pblh - REAL fak2(klon) ! k*wm*pblh - REAL fak3(klon) ! fakn*wstar/wm - REAL pblk(klon) ! level eddy diffusivity for momentum - REAL pr(klon) ! Prandtl number for eddy diffusivities - REAL zl(klon) ! zmzp / Obukhov length - REAL zh(klon) ! zmzp / pblh - REAL zzh(klon) ! (1-(zmzp/pblh))**2 - REAL zm(klon) ! current level height - REAL zp(klon) ! current level height + one level up - REAL zcor, zdelta, zcvm5 - ! AM REAL zxqs - REAL fac, pblmin, zmzp, term - - include "YOETHF.h" - include "FCTTRE.h" - - - - ! initialisations (Anne) - isommet = klev - th_th(:) = 0. - q_star = 0 - t_star = 0 - - - b212 = sqrt(b1*b2) - b2sr = sqrt(b2) - - ! ============================================================ - ! Fonctions thermo implicites - ! ============================================================ - ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - ! Tetens : pression partielle de vap d'eau e_sat(T) - ! ================================================= - ! ++ e_sat(T) = r2*exp( r3*(T-Tf)/(T-r4) ) id a r2*FOEWE - ! ++ avec : - ! ++ Tf = 273.16 K (Temp de fusion de la glace) - ! ++ r2 = 611.14 Pa - ! ++ r3 = 17.269 (liquide) 21.875 (solide) adim - ! ++ r4 = 35.86 7.66 Kelvin - ! ++ q_sat = eps*e_sat/(p-(1-eps)*e_sat) - ! ++ deriv� : - ! ++ ========= - ! ++ r3*(Tf-r4)*q_sat(T,p) - ! ++ d_qsat_dT = -------------------------------- - ! ++ (T-r4)^2*( 1-(1-eps)*e_sat(T)/p ) - ! ++ pour zcvm5=Lv, c'est FOEDE - ! ++ Rq :(1.-REPS)*esarg/Parg id a RETV*Qsat - ! ------------------------------------------------------------------ - - ! Initialisation - rlvcp = rlvtt/rcpd - reps = rd/rv - - - ! DO i = 1, klon - ! pcfh(i,1) = cd_h(i) - ! pcfm(i,1) = cd_m(i) - ! ENDDO - ! DO k = 2, klev - ! DO i = 1, klon - ! pcfh(i,k) = zkmin - ! pcfm(i,k) = zkmin - ! cgs(i,k) = 0.0 - ! cgh(i,k) = 0.0 - ! cgq(i,k) = 0.0 - ! ENDDO - ! ENDDO - - ! Calculer les hauteurs de chaque couche - ! (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g) - ! pourquoi ne pas utiliser Phi/RG ? - DO i = 1, knon - z(i, 1) = rd*t(i, 1)/(0.5*(paprs(i,1)+pplay(i,1)))*(paprs(i,1)-pplay(i,1) & - )/rg - s(i, 1) = (pplay(i,1)/paprs(i,1))**rkappa - END DO - ! s(k) = [pplay(k)/ps]^kappa - ! + + + + + + + + + pplay <-> s(k) t dp=pplay(k-1)-pplay(k) - - ! ----------------- paprs <-> sig(k) - - ! + + + + + + + + + pplay <-> s(k-1) - - - ! + + + + + + + + + pplay <-> s(1) t dp=paprs-pplay z(1) - - ! ----------------- paprs <-> sig(1) - - DO k = 2, klev - DO i = 1, knon - z(i, k) = z(i, k-1) + rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k)*(pplay(i,k-1 & - )-pplay(i,k))/rg - s(i, k) = (pplay(i,k)/paprs(i,1))**rkappa - END DO - END DO - ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - ! +++ Determination des grandeurs de surface +++++++++++++++++++++ - ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - DO i = 1, knon - ! AM IF (thermcep) THEN - ! AM zdelta=MAX(0.,SIGN(1.,RTT-tsol(i))) - ! zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta - ! zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1)) - ! AM zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1) - ! AM zxqs=MIN(0.5,zxqs) - ! AM zcor=1./(1.-retv*zxqs) - ! AM zxqs=zxqs*zcor - ! AM ELSE - ! AM IF (tsol(i).LT.t_coup) THEN - ! AM zxqs = qsats(tsol(i)) / paprs(i,1) - ! AM ELSE - ! AM zxqs = qsatl(tsol(i)) / paprs(i,1) - ! AM ENDIF - ! AM ENDIF - ! niveau de reference bulk; mais ici, c,a pourrait etre le niveau de ref - ! du thermique - ! AM zx_alf1 = 1.0 - ! AM zx_alf2 = 1.0 - zx_alf1 - ! AM zxt = (t(i,1)+z(i,1)*RG/RCPD/(1.+RVTMP2*q(i,1))) - ! AM . *(1.+RETV*q(i,1))*zx_alf1 - ! AM . + (t(i,2)+z(i,2)*RG/RCPD/(1.+RVTMP2*q(i,2))) - ! AM . *(1.+RETV*q(i,2))*zx_alf2 - ! AM zxu = u(i,1)*zx_alf1+u(i,2)*zx_alf2 - ! AM zxv = v(i,1)*zx_alf1+v(i,2)*zx_alf2 - ! AM zxq = q(i,1)*zx_alf1+q(i,2)*zx_alf2 - ! AM - ! AMAM zxu = u10m(i) - ! AMAM zxv = v10m(i) - ! AMAM zxmod = 1.0+SQRT(zxu**2+zxv**2) - ! AM Niveau de ref choisi a 2m - zxt = t2m(i) - - ! *************************************************** - ! attention, il doit s'agir de - ! ;Calcul de tcls virtuel et de w'theta'virtuel - ! ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; - ! tcls=tcls*(1+.608*qcls) - - ! ;Pour avoir w'theta', - ! ; il faut diviser par ro.Cp - ! Cp=Cpd*(1+0.84*qcls) - ! fcs=fcs/(ro_surf*Cp) - ! ;On transforme w'theta' en w'thetav' - ! Lv=(2.501-0.00237*(tcls-273.15))*1.E6 - ! xle=xle/(ro_surf*Lv) - ! fcsv=fcs+.608*xle*tcls - ! *************************************************** - ! AM khfs(i) = (tsol(i)*(1.+RETV*q(i,1))-zxt) *zxmod*cd_h(i) - ! AM kqfs(i) = (zxqs-zxq) *zxmod*cd_h(i) * beta(i) - ! AM - ! dif khfs est deja w't'_v / heatv(i) = khfs(i) + RETV*zxt*kqfs(i) - ! AM calcule de Ro = paprs(i,1)/Rd zxt - ! AM convention >0 vers le bas ds lmdz - khfs(i) = -flux_t(i, 1)*zxt*rd/(rcpd*paprs(i,1)) - kqfs(i) = -flux_q(i, 1)*zxt*rd/(paprs(i,1)) - ! AM verifier que khfs et kqfs sont bien de la forme w'l' - heatv(i) = khfs(i) + 0.608*zxt*kqfs(i) - ! a comparer aussi aux sorties de clqh : flux_T/RoCp et flux_q/RoLv - ! AM heatv(i) = khfs(i) - ! AM ustar est en entree - ! AM taux = zxu *zxmod*cd_m(i) - ! AM tauy = zxv *zxmod*cd_m(i) - ! AM ustar(i) = SQRT(taux**2+tauy**2) - ! AM ustar(i) = MAX(SQRT(ustar(i)),0.01) - ! Theta et qT du thermique sans exces (interpolin vers surf) - ! chgt de niveau du thermique (jeudi 30/12/1999) - ! (interpolation lineaire avant integration phi_h) - ! AM qT_th(i) = zxqs*beta(i) + 4./z(i,1)*(q(i,1)-zxqs*beta(i)) - ! AM qT_th(i) = max(qT_th(i),q(i,1)) - qt_th(i) = q2m(i) - ! n The_th restera la Theta du thermique sans exces jusqu'a 2eme calcul - ! n reste a regler convention P) pour Theta - ! The_th(i) = tsol(i) + 4./z(i,1)*(t(i,1)-tsol(i)) - ! - + RLvCp*qT_th(i) - ! AM Th_th(i) = tsol(i) + 4./z(i,1)*(t(i,1)-tsol(i)) - th_th(i) = t2m(i) - END DO - - DO i = 1, knon - rhino(i, 1) = 0.0 ! Global Richardson - check(i) = .TRUE. - pblh(i) = z(i, 1) ! on initialise pblh a l'altitude du 1er niveau - plcl(i) = 6000. - ! Lambda = -u*^3 / (alpha.g.kvon. - unsobklen(i) = -rg*vk*heatv(i)/(t(i,1)*max(ustar(i),usmin)**3) - trmb1(i) = 0. - trmb2(i) = 0. - trmb3(i) = 0. - END DO - - - ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - ! PBL height calculation: - ! Search for level of pbl. Scan upward until the Richardson number between - ! the first level and the current level exceeds the "critical" value. - ! (bonne idee Nu de separer le Ric et l'exces de temp du thermique) - ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - fac = 100.0 - DO k = 2, isommet - DO i = 1, knon - IF (check(i)) THEN - ! pourquoi / niveau 1 (au lieu du sol) et le terme en u*^2 ? - ! test zdu2 = - ! (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2 - zdu2 = u(i, k)**2 + v(i, k)**2 - zdu2 = max(zdu2, 1.0E-20) - ! Theta_v environnement - zthvd = t(i, k)/s(i, k)*(1.+retv*q(i,k)) - - ! therm Theta_v sans exces (avec hypothese fausse de H&B, sinon, - ! passer par Theta_e et virpot) - ! zthvu=t(i,1)/s(i,1)*(1.+RETV*q(i,1)) - ! AM zthvu = Th_th(i)*(1.+RETV*q(i,1)) - zthvu = th_th(i)*(1.+retv*qt_th(i)) - ! Le Ri par Theta_v - ! AM rhino(i,k) = (z(i,k)-z(i,1))*RG*(zthvd-zthvu) - ! AM . /(zdu2*0.5*(zthvd+zthvu)) - ! AM On a nveau de ref a 2m ??? - rhino(i, k) = (z(i,k)-zref)*rg*(zthvd-zthvu)/(zdu2*0.5*(zthvd+zthvu)) - - IF (rhino(i,k)>=ricr) THEN - pblh(i) = z(i, k-1) + (z(i,k-1)-z(i,k))*(ricr-rhino(i,k-1))/(rhino( & - i,k-1)-rhino(i,k)) - ! test04 - pblh(i) = pblh(i) + 100. - pblt(i) = t(i, k-1) + (t(i,k)-t(i,k-1))*(pblh(i)-z(i,k-1))/(z(i,k)- & - z(i,k-1)) - check(i) = .FALSE. - END IF - END IF - END DO - END DO - - - ! Set pbl height to maximum value where computation exceeds number of - ! layers allowed - - DO i = 1, knon - IF (check(i)) pblh(i) = z(i, isommet) - END DO - - ! Improve estimate of pbl height for the unstable points. - ! Find unstable points (sensible heat flux is upward): - - DO i = 1, knon - IF (heatv(i)>0.) THEN - unstbl(i) = .TRUE. - check(i) = .TRUE. - ELSE - unstbl(i) = .FALSE. - check(i) = .FALSE. - END IF - END DO - - ! For the unstable case, compute velocity scale and the - ! convective temperature excess: - - DO i = 1, knon - IF (check(i)) THEN - phiminv(i) = (1.-binm*pblh(i)*unsobklen(i))**onet - ! *************************************************** - ! Wm ? et W* ? c'est la formule pour z/h < .1 - ! ;Calcul de w* ;; - ! ;;;;;;;;;;;;;;;; - ! w_star=((g/tcls)*fcsv*z(ind))^(1/3.) [ou prendre la premiere approx - ! de h) - ! ;; CALCUL DE wm ;; - ! ;;;;;;;;;;;;;;;;;; - ! ; Ici on considerera que l'on est dans la couche de surf jusqu'a 100m - ! ; On prend svt couche de surface=0.1*h mais on ne connait pas h - ! ;;;;;;;;;;;Dans la couche de surface - ! if (z(ind) le 20) then begin - ! Phim=(1.-15.*(z(ind)/L))^(-1/3.) - ! wm=u_star/Phim - ! ;;;;;;;;;;;En dehors de la couche de surface - ! endif else if (z(ind) gt 20) then begin - ! wm=(u_star^3+c1*w_star^3)^(1/3.) - ! endif - ! *************************************************** - wm(i) = ustar(i)*phiminv(i) - ! ====================================================================== - ! valeurs de Dominique Lambert de la campagne SEMAPHORE : - ! = 100.T*^2; = 20.q*^2 a 10m - ! = (1+1.2q).100.T* + 1.2Tv.sqrt(20*100).T*.q* + - ! (.608*Tv)^2*20.q*^2; - ! et dTetavS = sqrt() ainsi calculee. - ! avec : T*=_s/w* et q*=/w* - ! !!! on peut donc utiliser w* pour les fluctuations <-> Lambert - ! (leur corellation pourrait dependre de beta par ex) - ! if fcsv(i,j) gt 0 then begin - ! dTetavs=b1*(1.+2.*.608*q_10(i,j))*(fcs(i,j)/wm(i,j))^2+$ - ! (.608*Thetav_10(i,j))^2*b2*(xle(i,j)/wm(i,j))^2+$ - ! 2.*.608*thetav_10(i,j)*sqrt(b1*b2)*(xle(i,j)/wm(i,j))*(fcs(i,j)/wm(i,j)) - ! dqs=b2*(xle(i,j)/wm(i,j))^2 - ! theta_s(i,j)=thetav_10(i,j)+sqrt(dTetavs) - ! q_s(i,j)=q_10(i,j)+sqrt(dqs) - ! endif else begin - ! Theta_s(i,j)=thetav_10(i,j) - ! q_s(i,j)=q_10(i,j) - ! endelse - ! ====================================================================== - - ! HBTM therm(i) = heatv(i)*fak/wm(i) - ! forme Mathieu : - q_star = kqfs(i)/wm(i) - t_star = khfs(i)/wm(i) - ! IM 091204 BEG - IF (1==0) THEN - IF (t_star<0. .OR. q_star<0.) THEN - PRINT *, 'i t_star q_star khfs kqfs wm', i, t_star, q_star, & - khfs(i), kqfs(i), wm(i) - END IF - END IF - ! IM 091204 END - ! AM Nveau cde ref 2m => - ! AM therm(i) = sqrt( b1*(1.+2.*RETV*q(i,1))*t_star**2 - ! AM + + (RETV*T(i,1))**2*b2*q_star**2 - ! AM + + 2.*RETV*T(i,1)*b212*q_star*t_star - ! AM + ) - ! IM 091204 BEG - a1 = b1*(1.+2.*retv*qt_th(i))*t_star**2 - a2 = (retv*th_th(i))**2*b2*q_star*q_star - a3 = 2.*retv*th_th(i)*b212*q_star*t_star - aa = a1 + a2 + a3 - IF (1==0) THEN - IF (aa<0.) THEN - PRINT *, 'i a1 a2 a3 aa', i, a1, a2, a3, aa - PRINT *, 'i qT_th Th_th t_star q_star RETV b1 b2 b212', i, & - qt_th(i), th_th(i), t_star, q_star, retv, b1, b2, b212 - END IF - END IF - ! IM 091204 END - therm(i) = sqrt(b1*(1.+2.*retv*qt_th(i))*t_star**2+(retv*th_th( & - i))**2*b2*q_star*q_star & ! IM 101204 + + - ! 2.*RETV*Th_th(i)*b212*q_star*t_star - +max(0.,2.*retv*th_th(i)*b212*q_star*t_star)) - - ! Theta et qT du thermique (forme H&B) avec exces - ! (attention, on ajoute therm(i) qui est virtuelle ...) - ! pourquoi pas sqrt(b1)*t_star ? - ! dqs = b2sr*kqfs(i)/wm(i) - qt_th(i) = qt_th(i) + b2sr*q_star - ! new on differre le calcul de Theta_e - ! The_th(i) = The_th(i) + therm(i) + RLvCp*qT_th(i) - ! ou: The_th(i) = The_th(i) + sqrt(b1)*khfs(i)/wm(i) + - ! RLvCp*qT_th(i) - rhino(i, 1) = 0.0 - END IF - END DO - - ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - ! ++ Improve pblh estimate for unstable conditions using the +++++++ - ! ++ convective temperature excess : +++++++ - ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - DO k = 2, isommet - DO i = 1, knon - IF (check(i)) THEN - ! test zdu2 = - ! (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2 - zdu2 = u(i, k)**2 + v(i, k)**2 - zdu2 = max(zdu2, 1.0E-20) - ! Theta_v environnement - zthvd = t(i, k)/s(i, k)*(1.+retv*q(i,k)) - - ! et therm Theta_v (avec hypothese de constance de H&B, - ! zthvu=(t(i,1)+therm(i))/s(i,1)*(1.+RETV*q(i,1)) - zthvu = th_th(i)*(1.+retv*qt_th(i)) + therm(i) - - - ! Le Ri par Theta_v - ! AM Niveau de ref 2m - ! AM rhino(i,k) = (z(i,k)-z(i,1))*RG*(zthvd-zthvu) - ! AM . /(zdu2*0.5*(zthvd+zthvu)) - rhino(i, k) = (z(i,k)-zref)*rg*(zthvd-zthvu)/(zdu2*0.5*(zthvd+zthvu)) - - - IF (rhino(i,k)>=ricr) THEN - pblh(i) = z(i, k-1) + (z(i,k-1)-z(i,k))*(ricr-rhino(i,k-1))/(rhino( & - i,k-1)-rhino(i,k)) - ! test04 - pblh(i) = pblh(i) + 100. - pblt(i) = t(i, k-1) + (t(i,k)-t(i,k-1))*(pblh(i)-z(i,k-1))/(z(i,k)- & - z(i,k-1)) - check(i) = .FALSE. - ! IM 170305 BEG - IF (1==0) THEN - ! debug print -120;34 -34- 58 et 0;26 wamp - IF (i==950 .OR. i==192 .OR. i==624 .OR. i==118) THEN - PRINT *, ' i,Th_th,Therm,qT :', i, th_th(i), therm(i), qt_th(i) - q_star = kqfs(i)/wm(i) - t_star = khfs(i)/wm(i) - PRINT *, 'q* t*, b1,b2,b212 ', q_star, t_star, & - b1*(1.+2.*retv*qt_th(i))*t_star**2, & - (retv*th_th(i))**2*b2*q_star**2, 2.*retv*th_th(i)*b212*q_star & - *t_star - PRINT *, 'zdu2 ,100.*ustar(i)**2', zdu2, fac*ustar(i)**2 - END IF - END IF !(1.EQ.0) THEN - ! IM 170305 END - ! q_star = kqfs(i)/wm(i) - ! t_star = khfs(i)/wm(i) - ! trmb1(i) = b1*(1.+2.*RETV*q(i,1))*t_star**2 - ! trmb2(i) = (RETV*T(i,1))**2*b2*q_star**2 - ! Omega now trmb3(i) = 2.*RETV*T(i,1)*b212*q_star*t_star - END IF - END IF - END DO - END DO - - ! Set pbl height to maximum value where computation exceeds number of - ! layers allowed - - DO i = 1, knon - IF (check(i)) pblh(i) = z(i, isommet) - END DO - - ! PBL height must be greater than some minimum mechanical mixing depth - ! Several investigators have proposed minimum mechanical mixing depth - ! relationships as a function of the local friction velocity, u*. We - ! make use of a linear relationship of the form h = c u* where c=700. - ! The scaling arguments that give rise to this relationship most often - ! represent the coefficient c as some constant over the local coriolis - ! parameter. Here we make use of the experimental results of Koracin - ! and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f - ! where f was evaluated at 39.5 N and 52 N. Thus we use a typical mid - ! latitude value for f so that c = 0.07/f = 700. - - DO i = 1, knon - pblmin = 700.0*ustar(i) - pblh(i) = max(pblh(i), pblmin) - ! par exemple : - pblt(i) = t(i, 2) + (t(i,3)-t(i,2))*(pblh(i)-z(i,2))/(z(i,3)-z(i,2)) - END DO - - ! ******************************************************************** - ! pblh is now available; do preparation for diffusivity calculation : - ! ******************************************************************** - DO i = 1, knon - check(i) = .TRUE. - zsat(i) = .FALSE. - ! omegafl utilise pour prolongement CAPE - omegafl(i) = .FALSE. - cape(i) = 0. - kape(i) = 0. - eauliq(i) = 0. - ctei(i) = 0. - pblk(i) = 0.0 - fak1(i) = ustar(i)*pblh(i)*vk - - ! Do additional preparation for unstable cases only, set temperature - ! and moisture perturbations depending on stability. - ! *** Rq: les formule sont prises dans leur forme CS *** - IF (unstbl(i)) THEN - ! AM Niveau de ref du thermique - ! AM zxt=(t(i,1)-z(i,1)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1))) - ! AM . *(1.+RETV*q(i,1)) - zxt = (th_th(i)-zref*0.5*rg/rcpd/(1.+rvtmp2*qt_th(i)))* & - (1.+retv*qt_th(i)) - phiminv(i) = (1.-binm*pblh(i)*unsobklen(i))**onet - phihinv(i) = sqrt(1.-binh*pblh(i)*unsobklen(i)) - wm(i) = ustar(i)*phiminv(i) - fak2(i) = wm(i)*pblh(i)*vk - wstar(i) = (heatv(i)*rg*pblh(i)/zxt)**onet - fak3(i) = fakn*wstar(i)/wm(i) - ELSE - wstar(i) = 0. - END IF - ! Computes Theta_e for thermal (all cases : to be modified) - ! attention ajout therm(i) = virtuelle - the_th(i) = th_th(i) + therm(i) + rlvcp*qt_th(i) - ! ou: The_th(i) = Th_th(i) + sqrt(b1)*khfs(i)/wm(i) + RLvCp*qT_th(i) - END DO - - ! Main level loop to compute the diffusivities and - ! counter-gradient terms: - - DO k = 2, isommet - - ! Find levels within boundary layer: - - DO i = 1, knon - unslev(i) = .FALSE. - stblev(i) = .FALSE. - zm(i) = z(i, k-1) - zp(i) = z(i, k) - IF (zkmin==0.0 .AND. zp(i)>pblh(i)) zp(i) = pblh(i) - IF (zm(i) 1e-3mm/day - - real ridicule_evap ! valeur limite de ridicule pour les evap - parameter (ridicule_evap=ridicule_rain*1e-2) ! en kg/s <-> 1e-3mm/day - - real ridicule_qsol ! valeur limite de ridicule pour les qsol - parameter (ridicule_qsol=ridicule_rain) ! en kg <-> 1e-8kg - - real ridicule_snow ! valeur limite de ridicule pour les snow - parameter (ridicule_snow=ridicule_qsol) ! en kg/s <-> 1e-8kg - - real expb_max - parameter (expb_max=30.0) - - ! spécifique au tritium: - - -logical, save :: ok_prod_nucl_tritium ! si oui, production de tritium par essais nucleaires -!$OMP THREADPRIVATE(ok_prod_nucl_tritium) - integer nessai - parameter (nessai=486) - integer, save :: day_nucl(nessai) -!$OMP THREADPRIVATE(day_nucl) - integer, save :: month_nucl(nessai) -!$OMP THREADPRIVATE(month_nucl) - integer, save :: year_nucl(nessai) -!$OMP THREADPRIVATE(year_nucl) - real, save :: lat_nucl(nessai) -!$OMP THREADPRIVATE(lat_nucl) - real, save :: lon_nucl(nessai) -!$OMP THREADPRIVATE(lon_nucl) - real, save :: zmin_nucl(nessai) -!$OMP THREADPRIVATE(zmin_nucl) - real, save :: zmax_nucl(nessai) -!$OMP THREADPRIVATE(zmax_nucl) - real, save :: HTO_nucl(nessai) -!$OMP THREADPRIVATE(HTO_nucl) - - -CONTAINS - - SUBROUTINE iso_init() - use IOIPSL ! getin - implicit none - -! -- local variables: - - integer ixt - ! référence O18 - real fac_enrichoce18 - real alpha_liq_sol_O18, & - & talph1_O18,talph2_O18,talph3_O18, & - & talps1_O18,talps2_O18, & - & tkcin0_O18,tkcin1_O18,tkcin2_O18, & - & tdifrel_O18 - - ! cas de l'O17 - real fac_kcin - real pente_MWL - integer ierr - - logical ok_nocinsat, ok_nocinsfc !sensi test - parameter (ok_nocinsfc=.FALSE.) ! if T: no kinetic effect in sfc evap - parameter (ok_nocinsat=.FALSE.) ! if T: no sursaturation effect for ice - logical Rdefault_smow - parameter (Rdefault_smow=.FALSE.) ! si T: Rdefault=smow; si F: nul - ! pour le tritium - integer iessai - -! allocations mémoire - -allocate (tnat(niso)) -allocate (toce(niso)) -allocate (tcorr(niso)) -allocate (tdifrel(niso)) -allocate (talph1(niso)) -allocate (talph2(niso)) -allocate (talph3(niso)) -allocate (talps1(niso)) -allocate (talps2(niso)) -allocate (tkcin0(niso)) -allocate (tkcin1(niso)) -allocate (tkcin2(niso)) -allocate (alpha_liq_sol(niso)) -allocate (Rdefault(niso)) -allocate (Rmethox(niso)) -allocate (striso(niso)) - - -!-------------------------------------------------------------- -! General: -!-------------------------------------------------------------- - -! -- verif du nombre d'isotopes: - write(*,*) 'iso_init 64: niso=',niso - -! init de ntracisoOR: on écrasera en cas de ok_isotrac si complications avec -! ORCHIDEE - ntracisoOR=ntraciso - -! -- Type of water isotopes: - - iso_eau=indnum_fn_num(1) - iso_HDO=indnum_fn_num(2) - iso_O18=indnum_fn_num(3) - iso_O17=indnum_fn_num(4) - iso_HTO=indnum_fn_num(5) - write(*,*) 'iso_init 59: iso_eau=',iso_eau - write(*,*) 'iso_HDO=',iso_HDO - write(*,*) 'iso_O18=',iso_O18 - write(*,*) 'iso_O17=',iso_O17 - write(*,*) 'iso_HTO=',iso_HTO - write(*,*) 'iso_init 251: use_iso=',use_iso - - ! initialisation - lambda_sursat=0.004 - thumxt1=0.75*1.2 - ntot=20 - h_land_ice=20. ! à comparer aux 3000mm de snow_max - P_veg=1.0 - bidouille_anti_divergence=.false. - essai_convergence=.false. - initialisation_iso=0 - modif_sst=0 - modif_sic=0 - deltaTtest=0.0 - deltasic=0.1 - deltaTtestpoles=0.0 - sstlatcrit=30.0 - deltaO18_oce=0.0 - albedo_prescrit=0 - lon_min_albedo=-200 - lon_max_albedo=200 - lat_min_albedo=-100 - lat_max_albedo=100 - deltaP_BL=10.0 - ruissellement_pluie=0 - alphak_stewart=1 - tdifexp_sol=0.67 - calendrier_guide=0 - cste_surf_cond=0 -mixlen=35.0 -evap_cont_cste=0.0 -deltaO18_evap_cont=0.0 -d_evap_cont=0.0 -nudge_qsol=0 -region_nudge_qsol=1 -nlevmaxO17=50 -no_pce=0 -A_satlim=1.0 -ok_restrict_A_satlim=0 -! slope_limiterxy=2.0 -! slope_limiterz=2.0 -modif_ratqs=0 -Pcrit_ratqs=500.0 -ratqsbasnew=0.05 - -fac_modif_evaoce=1.0 -ok_bidouille_wake=0 -cond_temp_env=.false. -! si oui, la temperature de cond est celle de l'environnement, -! pour eviter bugs quand temperature dans ascendances convs est -! mal calculee -ok_prod_nucl_tritium=.false. - -! lecture des paramètres isotopiques: -call getin('lambda',lambda_sursat) -call getin('thumxt1',thumxt1) -call getin('ntot',ntot) -call getin('h_land_ice',h_land_ice) -call getin('P_veg',P_veg) -call getin('bidouille_anti_divergence',bidouille_anti_divergence) -call getin('essai_convergence',essai_convergence) -call getin('initialisation_iso',initialisation_iso) -!if (ok_isotrac) then -!if (initialisation_iso.eq.0) then -! call getin('initialisation_isotrac',initialisation_isotrac) -!endif !if (initialisation_iso.eq.0) then -!endif !if (ok_isotrac) then -call getin('modif_sst',modif_sst) -if (modif_sst.ge.1) then -call getin('deltaTtest',deltaTtest) -if (modif_sst.ge.2) then - call getin('deltaTtestpoles',deltaTtestpoles) - call getin('sstlatcrit',sstlatcrit) -#ifdef ISOVERIF - !call iso_verif_positif(sstlatcrit,'iso_init 107') - if (sstlatcrit.lt.0.0) then - write(*,*) 'iso_init 270: sstlatcrit=',sstlatcrit - stop - endif -#endif - if (modif_sst.ge.3) then - call getin('dsstlatcrit',dsstlatcrit) -#ifdef ISOVERIF - !call iso_verif_positif(dsstlatcrit,'iso_init 110') - if (sstlatcrit.lt.0.0) then - write(*,*) 'iso_init 279: dsstlatcrit=',dsstlatcrit - stop - endif -#endif - endif !if (modif_sst.ge.3) then -endif !if (modif_sst.ge.2) then -endif ! if (modif_sst.ge.1) then -call getin('modif_sic',modif_sic) -if (modif_sic.ge.1) then -call getin('deltasic',deltasic) -endif !if (modif_sic.ge.1) then - -call getin('albedo_prescrit',albedo_prescrit) -call getin('lon_min_albedo',lon_min_albedo) -call getin('lon_max_albedo',lon_max_albedo) -call getin('lat_min_albedo',lat_min_albedo) -call getin('lat_max_albedo',lat_max_albedo) -call getin('deltaO18_oce',deltaO18_oce) -call getin('deltaP_BL',deltaP_BL) -call getin('ruissellement_pluie',ruissellement_pluie) -call getin('alphak_stewart',alphak_stewart) -call getin('tdifexp_sol',tdifexp_sol) -call getin('calendrier_guide',calendrier_guide) -call getin('cste_surf_cond',cste_surf_cond) -call getin('mixlen',mixlen) -call getin('evap_cont_cste',evap_cont_cste) -call getin('deltaO18_evap_cont',deltaO18_evap_cont) -call getin('d_evap_cont',d_evap_cont) -call getin('nudge_qsol',nudge_qsol) -call getin('region_nudge_qsol',region_nudge_qsol) -call getin('no_pce',no_pce) -call getin('A_satlim',A_satlim) -call getin('ok_restrict_A_satlim',ok_restrict_A_satlim) -#ifdef ISOVERIF -!call iso_verif_positif(1.0-A_satlim,'iso_init 158') - if (A_satlim.gt.1.0) then - write(*,*) 'iso_init 315: A_satlim=',A_satlim - stop - endif -#endif -! call getin('slope_limiterxy',slope_limiterxy) -! call getin('slope_limiterz',slope_limiterz) -call getin('modif_ratqs',modif_ratqs) -call getin('Pcrit_ratqs',Pcrit_ratqs) -call getin('ratqsbasnew',ratqsbasnew) -call getin('fac_modif_evaoce',fac_modif_evaoce) -call getin('ok_bidouille_wake',ok_bidouille_wake) -call getin('cond_temp_env',cond_temp_env) -if (use_iso(iso_HTO_possible)) then - ok_prod_nucl_tritium=.true. - call getin('ok_prod_nucl_tritium',ok_prod_nucl_tritium) -endif - -write(*,*) 'lambda,thumxt1=',lambda_sursat,thumxt1 -write(*,*) 'bidouille_anti_divergence=',bidouille_anti_divergence -write(*,*) 'essai_convergence=',essai_convergence -write(*,*) 'initialisation_iso=',initialisation_iso -write(*,*) 'modif_sst=',modif_sst -if (modif_sst.ge.1) then -write(*,*) 'deltaTtest=',deltaTtest -if (modif_sst.ge.2) then -write(*,*) 'deltaTtestpoles,sstlatcrit=', & -& deltaTtestpoles,sstlatcrit -if (modif_sst.ge.3) then - write(*,*) 'dsstlatcrit=',dsstlatcrit -endif !if (modif_sst.ge.3) then -endif !if (modif_sst.ge.2) then -endif !if (modif_sst.ge.1) then -write(*,*) 'modif_sic=',modif_sic -if (modif_sic.ge.1) then -write(*,*) 'deltasic=',deltasic -endif !if (modif_sic.ge.1) then -write(*,*) 'deltaO18_oce=',deltaO18_oce -write(*,*) 'albedo_prescrit=',albedo_prescrit -if (albedo_prescrit.eq.1) then - write(*,*) 'lon_min_albedo,lon_max_albedo=', & -& lon_min_albedo,lon_max_albedo - write(*,*) 'lat_min_albedo,lat_max_albedo=', & -& lat_min_albedo,lat_max_albedo -endif !if (albedo_prescrit.eq.1) then -write(*,*) 'deltaP_BL,ruissellement_pluie,alphak_stewart=', & -& deltaP_BL,ruissellement_pluie,alphak_stewart -write(*,*) 'cste_surf_cond=',cste_surf_cond -write(*,*) 'mixlen=',mixlen -write(*,*) 'tdifexp_sol=',tdifexp_sol -write(*,*) 'calendrier_guide=',calendrier_guide -write(*,*) 'evap_cont_cste=',evap_cont_cste -write(*,*) 'deltaO18_evap_cont,d_evap_cont=', & -& deltaO18_evap_cont,d_evap_cont -write(*,*) 'nudge_qsol,region_nudge_qsol=', & -& nudge_qsol,region_nudge_qsol -write(*,*) 'nlevmaxO17=',nlevmaxO17 -write(*,*) 'no_pce=',no_pce -write(*,*) 'A_satlim=',A_satlim -write(*,*) 'ok_restrict_A_satlim=',ok_restrict_A_satlim -! write(*,*) 'slope_limiterxy=',slope_limiterxy -! write(*,*) 'slope_limiterz=',slope_limiterz -write(*,*) 'modif_ratqs=',modif_ratqs -write(*,*) 'Pcrit_ratqs=',Pcrit_ratqs -write(*,*) 'ratqsbasnew=',ratqsbasnew -write(*,*) 'fac_modif_evaoce=',fac_modif_evaoce -write(*,*) 'ok_bidouille_wake=',ok_bidouille_wake -write(*,*) 'cond_temp_env=',cond_temp_env -write(*,*) 'ok_prod_nucl_tritium=',ok_prod_nucl_tritium - - -!-------------------------------------------------------------- -! Parameters that do not depend on the nature of water isotopes: -!-------------------------------------------------------------- - -! -- temperature at which ice condensate starts to form (valeur ECHAM?): -pxtmelt=273.15 -! pxtmelt=273.15-10.0 ! test PHASE - -! -- temperature at which all condensate is ice: -pxtice=273.15-10.0 -! pxtice=273.15-30.0 ! test PHASE - -! -- minimum temperature to calculate fractionation coeff -pxtmin=273.15-120.0 ! On ne calcule qu'au dessus de -120°C -pxtmax=273.15+60.0 ! On ne calcule qu'au dessus de +60°C -! remarque: les coeffs ont été mesurés seulement jusq'à -40! - -! -- a constant for alpha_eff for equilibrium below cloud base: -tdifexp=0.58 -tv0cin=7.0 - -! facteurs lambda et mu dans Si=musi-lambda*T -musi=1.0 -if (ok_nocinsat) then -lambda_sursat = 0.0 ! no sursaturation effect -endif - - -! diffusion dans le sol -Kd=2.5e-9 ! m2/s - -! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir -rh_cste_surf_cond=0.6 -T_cste_surf_cond=288.0 - -!-------------------------------------------------------------- -! Parameters that depend on the nature of water isotopes: -!-------------------------------------------------------------- -! ** constantes locales -fac_enrichoce18=0.0005 -! on a alors tcor018=1+fac_enrichoce18 -! tcorD=1+fac_enrichoce18*8 -! tcorO17=1+fac_enrichoce18*0.528 -alpha_liq_sol_O18=1.00291 ! valeur de Lehmann & Siegenthaler, 1991, - ! Journal of Glaciology, vol 37, p 23 -talph1_O18=1137. -talph2_O18=-0.4156 -talph3_O18=-2.0667E-3 -talps1_O18=11.839 -talps2_O18=-0.028244 -tkcin0_O18 = 0.006 -tkcin1_O18 = 0.000285 -tkcin2_O18 = 0.00082 -tdifrel_O18= 1./0.9723 - -! rapport des ln(alphaeq) entre O18 et O17 -fac_coeff_eq17_liq=0.529 ! donné par Amaelle -! fac_coeff_eq17_ice=0.528 ! slope MWL -fac_coeff_eq17_ice=0.529 - - -write(*,*) 'iso_O18,iso_HDO,iso_eau=',iso_O18,iso_HDO,iso_eau -do 999 ixt = 1, niso -write(*,*) 'iso_init 80: ixt=',ixt - - -! -- kinetic factor for surface evaporation: -! (cf: kcin = tkcin0 if |V|tv0cin ) -! (Rq: formula discontinuous for |V|=tv0cin... ) - -! -- main: -if (ixt.eq.iso_HTO) then ! Tritium - tkcin0(ixt) = 0.01056 - tkcin1(ixt) = 0.0005016 - tkcin2(ixt) = 0.0014432 - tnat(ixt)=0. - !toce(ixt)=2.2222E-8 ! corrigé par Alex Cauquoin - !toce(ixt)=1.0E-18 ! rapport 3H/1H ocean - toce(ixt)=4.0E-19 ! rapport T/H = 0.2 TU Dreisigacker and Roether 1978 - tcorr(ixt)=1. - tdifrel(ixt)=1./0.968 - talph1(ixt)=46480. - talph2(ixt)=-103.87 - talph3(ixt)=0. - talps1(ixt)=46480. - talps2(ixt)=-103.87 - alpha_liq_sol(ixt)=1. - Rdefault(ixt)=0.0 - Rmethox(ixt)=0.0 - striso(ixt)='HTO' -endif -if (ixt.eq.iso_O17) then ! Deuterium - pente_MWL=0.528 -! tdifrel(ixt)=1./0.985452 ! donné par Amaelle - tdifrel(ixt)=1./0.98555 ! valeur utilisée en 1D et dans modèle de LdG -! fac_kcin=0.5145 ! donné par Amaelle - fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0) - tkcin0(ixt) = tkcin0_O18*fac_kcin - tkcin1(ixt) = tkcin1_O18*fac_kcin - tkcin2(ixt) = tkcin2_O18*fac_kcin - tnat(ixt)=0.004/100. ! O17 représente 0.004% de l'oxygène - toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)**pente_MWL - tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL ! donné par Amaelle - talph1(ixt)=talph1_O18 - talph2(ixt)=talph2_O18 - talph3(ixt)=talph3_O18 - talps1(ixt)=talps1_O18 - talps2(ixt)=talps2_O18 - alpha_liq_sol(ixt)=(alpha_liq_sol_O18)**fac_coeff_eq17_liq - if (Rdefault_smow) then - Rdefault(ixt)=tnat(ixt)*(-3.15/1000.0+1.0) - else - Rdefault(ixt)=0.0 - endif - Rmethox(ixt)=(230./1000.+1.)*tnat(ixt) !Zahn et al 2006 - striso(ixt)='O17' -endif - -if (ixt.eq.iso_O18) then ! Oxygene18 - tkcin0(ixt) = tkcin0_O18 - tkcin1(ixt) = tkcin1_O18 - tkcin2(ixt) = tkcin2_O18 - tnat(ixt)=2005.2E-6 - toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0) - tcorr(ixt)=1.0+fac_enrichoce18 - tdifrel(ixt)=tdifrel_O18 - talph1(ixt)=talph1_O18 - talph2(ixt)=talph2_O18 - talph3(ixt)=talph3_O18 - talps1(ixt)=talps1_O18 - talps2(ixt)=talps2_O18 - alpha_liq_sol(ixt)=alpha_liq_sol_O18 - if (Rdefault_smow) then - Rdefault(ixt)=tnat(ixt)*(-6.0/1000.0+1.0) - else - Rdefault(ixt)=0.0 - endif - Rmethox(ixt)=(130./1000.+1.)*tnat(ixt) !Zahn et al 2006 -! write(*,*) 'iso_init 163: ZXalpha_liq_sol=',ZXalpha_liq_sol - striso(ixt)='O18' - write(*,*) 'isotopes_mod 519: ixt,striso(ixt)=',ixt,striso(ixt) -endif - -if (ixt.eq.iso_HDO) then ! Deuterium - pente_MWL=8.0 -! fac_kcin=0.88 - tdifrel(ixt)=1./0.9755 - fac_kcin= (tdifrel(ixt)-1)/(tdifrel_O18-1) - tkcin0(ixt) = tkcin0_O18*fac_kcin - tkcin1(ixt) = tkcin1_O18*fac_kcin - tkcin2(ixt) = tkcin2_O18*fac_kcin - tnat(ixt)=155.76E-6 - toce(ixt)=tnat(ixt)*(1.0+pente_MWL*deltaO18_oce/1000.0) - tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL - talph1(ixt)=24844. - talph2(ixt)=-76.248 - talph3(ixt)=52.612E-3 - talps1(ixt)=16288. - talps2(ixt)=-0.0934 - !ZXalpha_liq_sol=1.0192 ! Weston, Ralph, 1955 - alpha_liq_sol(ixt)=1.0212 - ! valeur de Lehmann & Siegenthaler, 1991, Journal of - ! Glaciology, vol 37, p 23 - if (Rdefault_smow) then - Rdefault(ixt)=tnat(ixt)*((-6.0*pente_MWL+10.0)/1000.0+1.0) - else - Rdefault(ixt)=0.0 - endif -Rmethox(ixt)=tnat(ixt)*(-25.0/1000.+1.) ! Zahn et al 2006 -striso(ixt)='HDO' - write(*,*) 'isotopes_mod 548: ixt,striso(ixt)=',ixt,striso(ixt) -endif - -! write(*,*) 'iso_init 163: ZXalpha_liq_sol=',ZXalpha_liq_sol -if (ixt.eq.iso_eau) then ! Oxygene16 - tkcin0(ixt) = 0.0 - tkcin1(ixt) = 0.0 - tkcin2(ixt) = 0.0 - tnat(ixt)=1. - toce(ixt)=tnat(ixt) - tcorr(ixt)=1.0 - tdifrel(ixt)=1. - talph1(ixt)=0. - talph2(ixt)=0. - talph3(ixt)=0. - talps1(ixt)=0. - talph3(ixt)=0. - alpha_liq_sol(ixt)=1. - if (Rdefault_smow) then - Rdefault(ixt)=tnat(ixt)*1.0 - else - Rdefault(ixt)=1.0 - endif - Rmethox(ixt)=1.0 - striso(ixt)='eau' -endif - -999 continue - -! test de sensibilité: -if (ok_nocinsfc) then ! no kinetic effect in sfc evaporation - do ixt=1,niso - tkcin0(ixt) = 0.0 - tkcin1(ixt) = 0.0 - tkcin2(ixt) = 0.0 - enddo -endif - -! fermeture fichier de paramètres -close(unit=32) - -! nom des isotopes - -! verif -write(*,*) 'iso_init 285: verif initialisation:' - -do ixt=1,niso - write(*,*) '* striso(',ixt,')=<'//striso(ixt)//'>' - write(*,*) 'tnat(',ixt,')=',tnat(ixt) -! write(*,*) 'alpha_liq_sol(',ixt,')=',alpha_liq_sol(ixt) -! write(*,*) 'tkcin0(',ixt,')=',tkcin0(ixt) -! write(*,*) 'tdifrel(',ixt,')=',tdifrel(ixt) -enddo -write(*,*) 'iso_init 69: lambda=',lambda_sursat -write(*,*) 'iso_init 69: thumxt1=',thumxt1 -write(*,*) 'iso_init 69: h_land_ice=',h_land_ice -write(*,*) 'iso_init 69: P_veg=',P_veg - - -END SUBROUTINE iso_init - -! functions basiques mises ici pour éviter dépendances circulaires - !*********** - function double_to_real(a) - implicit none - double precision a - real double_to_real - - double_to_real=a - - return - end function double_to_real - - !*********** - function real_to_double(a) - implicit none - real a - double precision real_to_double - - real_to_double=a - - return - end function real_to_double - -END MODULE isotopes_mod -#endif - - Index: LMDZ6/trunk/libf/phylmdiso/screenc.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/screenc.F90 (revision 3940) +++ (revision ) @@ -1,96 +1,0 @@ -! -! $Header$ -! - SUBROUTINE screenc(klon, knon, nsrf, zxli, & - speed, temp, q_zref, zref, & - ts, qsurf, z0m, z0h, psol, & - ustar, testar, qstar, okri, ri1, & - pref, delu, delte, delq) - IMPLICIT NONE -!----------------------------------------------------------------------- -! -! Objet : calcul "correcteur" des anomalies du vent, de la temperature -! potentielle et de l'humidite relative au niveau de reference zref et -! par rapport au 1er niveau (pour u) ou a la surface (pour theta et q) -! a partir des equations de Louis. -! -! Reference : Hess, Colman et McAvaney (1995) -! -! I. Musat, 01.07.2002 -!----------------------------------------------------------------------- -! -! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude) -! knon----input-I- nombre de points pour un type de surface -! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90 -! zxli----input-L- TRUE si calcul des cdrags selon Laurent Li -! speed---input-R- module du vent au 1er niveau du modele -! temp----input-R- temperature de l'air au 1er niveau du modele -! q_zref--input-R- humidite relative au 1er niveau du modele -! zref----input-R- altitude de reference -! ts------input-R- temperature de l'air a la surface -! qsurf---input-R- humidite relative a la surface -! z0m, z0h---input-R- rugosite -! psol----input-R- pression au sol -! ustar---input-R- facteur d'echelle pour le vent -! testar--input-R- facteur d'echelle pour la temperature potentielle -! qstar---input-R- facteur d'echelle pour l'humidite relative -! okri----input-L- TRUE si on veut tester le nb. Richardson entre la sfce -! et zref par rapport au Ri entre la sfce et la 1ere couche -! ri1-----input-R- nb. Richardson entre la surface et la 1ere couche -! -! pref----input-R- pression au niveau de reference -! delu----input-R- anomalie du vent par rapport au 1er niveau -! delte---input-R- anomalie de la temperature potentielle par rapport a la surface -! delq----input-R- anomalie de l'humidite relative par rapport a la surface -! - INTEGER, intent(in) :: klon, knon, nsrf - LOGICAL, intent(in) :: zxli, okri - REAL, dimension(klon), intent(in) :: speed, temp, q_zref - REAL, intent(in) :: zref - REAL, dimension(klon), intent(in) :: ts, qsurf, z0m, z0h, psol - REAL, dimension(klon), intent(in) :: ustar, testar, qstar, ri1 -! - REAL, dimension(klon), intent(out) :: pref, delu, delte, delq -!----------------------------------------------------------------------- - include "YOMCST.h" - include "flux_arp.h" -! -! Variables locales - INTEGER :: i - REAL, dimension(klon) :: cdram, cdrah, cdran, zri1, gref,ycdragm -! -!------------------------------------------------------------------------- - DO i=1, knon - gref(i) = zref*RG - ENDDO -! -! Richardson at reference level -! -! CALL coefcdrag (klon, knon, nsrf, zxli, & -! speed, temp, q_zref, gref, & -! psol, ts, qsurf, rugos, & -! okri, ri1, & -! cdram, cdrah, cdran, zri1, & -! pref) -! Fuxing WANG, 04/03/2015, replace the coefcdrag by the merged version: cdrag - CALL cdrag (knon, nsrf, & - speed, temp, q_zref, gref, & - psol, ts, qsurf, z0m, z0h, & - cdram, cdrah, zri1, pref) - DO i = 1, knon - IF(ok_prescr_ust) THEN -! La aussi il faut forcer avec ust (FC + MPL 20160210) - ycdragm(i) = ust*ust/(1.+speed(i))/speed(i) - cdram=ycdragm - delu(i) = ust/sqrt(cdram(i)) - ELSE - delu(i) = ustar(i)/sqrt(cdram(i)) - ENDIF - delte(i)= (testar(i)* sqrt(cdram(i)))/ & - cdrah(i) - delq(i)= (qstar(i)* sqrt(cdram(i)))/ & - cdrah(i) - ENDDO -! - RETURN - END SUBROUTINE screenc Index: LMDZ6/trunk/libf/phylmdiso/screenp.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/screenp.F90 (revision 3940) +++ (revision ) @@ -1,109 +1,0 @@ -! -! $Header$ -! - SUBROUTINE screenp(klon, knon, nsrf, & - & speed, tair, qair, & - & ts, qsurf, rugos, lmon, & - & ustar, testar, qstar, zref, & - & delu, delte, delq) - IMPLICIT none -!------------------------------------------------------------------------- -! -! Objet : calcul "predicteur" des anomalies du vent, de la temperature -! potentielle et de l'humidite relative au niveau de reference zref et -! par rapport au 1er niveau (pour u) ou a la surface (pour theta et q) -! a partir des relations de Dyer-Businger. -! -! Reference : Hess, Colman et McAvaney (1995) -! -! I. Musat, 01.07.2002 -!------------------------------------------------------------------------- -! -! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude) -! knon----input-I- nombre de points pour un type de surface -! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90 -! speed---input-R- module du vent au 1er niveau du modele -! tair----input-R- temperature de l'air au 1er niveau du modele -! qair----input-R- humidite relative au 1er niveau du modele -! ts------input-R- temperature de l'air a la surface -! qsurf---input-R- humidite relative a la surface -! rugos---input-R- rugosite -! lmon----input-R- longueur de Monin-Obukov -! ustar---input-R- facteur d'echelle pour le vent -! testar--input-R- facteur d'echelle pour la temperature potentielle -! qstar---input-R- facteur d'echelle pour l'humidite relative -! zref----input-R- altitude de reference -! -! delu----input-R- anomalie du vent par rapport au 1er niveau -! delte---input-R- anomalie de la temperature potentielle par rapport a la surface -! delq----input-R- anomalie de l'humidite relative par rapport a la surface -! - INTEGER, intent(in) :: klon, knon, nsrf - REAL, dimension(klon), intent(in) :: speed, tair, qair - REAL, dimension(klon), intent(in) :: ts, qsurf, rugos - DOUBLE PRECISION, dimension(klon), intent(in) :: lmon - REAL, dimension(klon), intent(in) :: ustar, testar, qstar - REAL, intent(in) :: zref -! - REAL, dimension(klon), intent(out) :: delu, delte, delq -! -!------------------------------------------------------------------------- -! Variables locales et constantes : - REAL, PARAMETER :: RKAR=0.40 - INTEGER :: i - REAL :: xtmp, xtmp0 -!------------------------------------------------------------------------- - DO i = 1, knon -! - IF (lmon(i).GE.0.) THEN -! -! STABLE CASE -! - IF (speed(i).GT.1.5.AND.lmon(i).LE.1.0 & - & .AND. rugos(i).LE.1.0) THEN - delu(i) = (ustar(i)/RKAR)* & - (log(zref/(rugos(i))+1.) + & - min(5.d0, 5.0 *(zref - rugos(i))/lmon(i))) - delte(i) = (testar(i)/RKAR)* & - (log(zref/(rugos(i))+1.) + & - min(5.d0, 5.0 * (zref - rugos(i))/lmon(i))) - delq(i) = (qstar(i)/RKAR)* & - (log(zref/(rugos(i))+1.) + & - min(5.d0, 5.0 * (zref - rugos(i))/lmon(i))) - ELSE - delu(i) = 0.1 * speed(i) - delte(i) = 0.1 * (tair(i) - ts(i) ) - delq(i) = 0.1 * (max(qair(i),0.0) - max(qsurf(i),0.0)) - ENDIF - ELSE -! -! UNSTABLE CASE -! - IF (speed(i).GT.5.0.AND.abs(lmon(i)).LE.50.0) THEN - xtmp = (1. - 16. * (zref/lmon(i)))**(1./4.) - xtmp0 = (1. - 16. * (rugos(i)/lmon(i)))**(1./4.) - delu(i) = (ustar(i)/RKAR)* & - (log(zref/(rugos(i))+1.) & - - 2.*log(0.5*(1. + xtmp)) & - + 2.*log(0.5*(1. + xtmp0)) & - - log(0.5*(1. + xtmp*xtmp)) & - + log(0.5*(1. + xtmp0*xtmp0)) & - + 2.*atan(xtmp) - 2.*atan(xtmp0)) - delte(i) = (testar(i)/RKAR)* & - (log(zref/(rugos(i))+1.) & - - 2.0 * log(0.5*(1. + xtmp*xtmp)) & - + 2.0 * log(0.5*(1. + xtmp0*xtmp0))) - delq(i) = (qstar(i)/RKAR)* & - (log(zref/(rugos(i))+1.) & - - 2.0 * log(0.5*(1. + xtmp*xtmp)) & - + 2.0 * log(0.5*(1. + xtmp0*xtmp0))) - ELSE - delu(i) = 0.5 * speed(i) - delte(i) = 0.5 * (tair(i) - ts(i) ) - delq(i) = 0.5 * (max(qair(i),0.0) - max(qsurf(i),0.0)) - ENDIF - ENDIF -! - ENDDO - RETURN - END SUBROUTINE screenp Index: LMDZ6/trunk/libf/phylmdiso/stdlevvar.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/stdlevvar.F90 (revision 3940) +++ (revision ) @@ -1,294 +1,0 @@ -! -! $Header$ -! - SUBROUTINE stdlevvar(klon, knon, nsrf, zxli, & - u1, v1, t1, q1, z1, & - ts1, qsurf, z0m, z0h, psol, pat1, & - t_2m, q_2m, t_10m, q_10m, u_10m, ustar) - IMPLICIT NONE -!------------------------------------------------------------------------- -! -! Objet : calcul de la temperature et l'humidite relative a 2m et du -! module du vent a 10m a partir des relations de Dyer-Businger et -! des equations de Louis. -! -! Reference : Hess, Colman et McAvaney (1995) -! -! I. Musat, 01.07.2002 -! -!AM On rajoute en sortie t et q a 10m pr le calcule d'hbtm2 dans clmain -! -!------------------------------------------------------------------------- -! -! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude) -! knon----input-I- nombre de points pour un type de surface -! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90 -! zxli----input-L- TRUE si calcul des cdrags selon Laurent Li -! u1------input-R- vent zonal au 1er niveau du modele -! v1------input-R- vent meridien au 1er niveau du modele -! t1------input-R- temperature de l'air au 1er niveau du modele -! q1------input-R- humidite relative au 1er niveau du modele -! z1------input-R- geopotentiel au 1er niveau du modele -! ts1-----input-R- temperature de l'air a la surface -! qsurf---input-R- humidite relative a la surface -! z0m, z0h---input-R- rugosite -! psol----input-R- pression au sol -! pat1----input-R- pression au 1er niveau du modele -! -! t_2m---output-R- temperature de l'air a 2m -! q_2m---output-R- humidite relative a 2m -! u_10m--output-R- vitesse du vent a 10m -!AM -! t_10m--output-R- temperature de l'air a 10m -! q_10m--output-R- humidite specifique a 10m -! ustar--output-R- u* -! - INTEGER, intent(in) :: klon, knon, nsrf - LOGICAL, intent(in) :: zxli - REAL, dimension(klon), intent(in) :: u1, v1, t1, q1, z1, ts1 - REAL, dimension(klon), intent(in) :: qsurf, z0m, z0h - REAL, dimension(klon), intent(in) :: psol, pat1 -! - REAL, dimension(klon), intent(out) :: t_2m, q_2m, ustar - REAL, dimension(klon), intent(out) :: u_10m, t_10m, q_10m -!------------------------------------------------------------------------- - include "flux_arp.h" - include "YOMCST.h" -!IM PLUS - include "YOETHF.h" -! -! Quelques constantes et options: -! -! RKAR : constante de von Karman - REAL, PARAMETER :: RKAR=0.40 -! niter : nombre iterations calcul "corrector" -! INTEGER, parameter :: niter=6, ncon=niter-1 - INTEGER, parameter :: niter=2, ncon=niter-1 -! -! Variables locales - INTEGER :: i, n - REAL :: zref - REAL, dimension(klon) :: speed -! tpot : temperature potentielle - REAL, dimension(klon) :: tpot - REAL, dimension(klon) :: zri1, cdran - REAL, dimension(klon) :: cdram, cdrah -! ri1 : nb. de Richardson entre la surface --> la 1ere couche - REAL, dimension(klon) :: ri1 - REAL, dimension(klon) :: testar, qstar - REAL, dimension(klon) :: zdte, zdq -! lmon : longueur de Monin-Obukhov selon Hess, Colman and McAvaney - DOUBLE PRECISION, dimension(klon) :: lmon - DOUBLE PRECISION, parameter :: eps=1.0D-20 - REAL, dimension(klon) :: delu, delte, delq - REAL, dimension(klon) :: u_zref, te_zref, q_zref - REAL, dimension(klon) :: temp, pref - LOGICAL :: okri - REAL, dimension(klon) :: u_zref_p, te_zref_p, temp_p, q_zref_p -!convertgence - REAL, dimension(klon) :: te_zref_con, q_zref_con - REAL, dimension(klon) :: u_zref_c, te_zref_c, temp_c, q_zref_c - REAL, dimension(klon) :: ok_pred, ok_corr -! REAL, dimension(klon) :: conv_te, conv_q -!------------------------------------------------------------------------- - DO i=1, knon - speed(i)=SQRT(u1(i)**2+v1(i)**2) - ri1(i) = 0.0 - ENDDO -! - okri=.FALSE. -! CALL coefcdrag(klon, knon, nsrf, zxli, & -! & speed, t1, q1, z1, psol, & -! & ts1, qsurf, rugos, okri, ri1, & -! & cdram, cdrah, cdran, zri1, pref) -! Fuxing WANG, 04/03/2015, replace the coefcdrag by the merged version: cdrag - - CALL cdrag(knon, nsrf, & - & speed, t1, q1, z1, & - & psol, ts1, qsurf, z0m, z0h, & - & cdram, cdrah, zri1, pref) - -! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013 - IF (ok_prescr_ust) then - DO i = 1, knon - print *,'cdram avant=',cdram(i) - cdram(i) = ust*ust/speed(i)/speed(i) - print *,'cdram ust speed apres=',cdram(i),ust,speed - ENDDO - ENDIF -! -!---------Star variables---------------------------------------------------- -! - DO i = 1, knon - ri1(i) = zri1(i) - tpot(i) = t1(i)* (psol(i)/pat1(i))**RKAPPA - ustar(i) = sqrt(cdram(i) * speed(i) * speed(i)) - zdte(i) = tpot(i) - ts1(i) - zdq(i) = max(q1(i),0.0) - max(qsurf(i),0.0) -! -! -!IM BUG BUG BUG zdte(i) = max(zdte(i),1.e-10) - zdte(i) = sign(max(abs(zdte(i)),1.e-10),zdte(i)) -! - testar(i) = (cdrah(i) * zdte(i) * speed(i))/ustar(i) - qstar(i) = (cdrah(i) * zdq(i) * speed(i))/ustar(i) - lmon(i) = (ustar(i) * ustar(i) * tpot(i))/ & - & (RKAR * RG * testar(i)) - ENDDO -! -!----------First aproximation of variables at zref -------------------------- - zref = 2.0 - CALL screenp(klon, knon, nsrf, speed, tpot, q1, & - & ts1, qsurf, z0m, lmon, & - & ustar, testar, qstar, zref, & - & delu, delte, delq) -! - DO i = 1, knon - u_zref(i) = delu(i) - q_zref(i) = max(qsurf(i),0.0) + delq(i) - te_zref(i) = ts1(i) + delte(i) - temp(i) = te_zref(i) * (psol(i)/pat1(i))**(-RKAPPA) - q_zref_p(i) = q_zref(i) -! te_zref_p(i) = te_zref(i) - temp_p(i) = temp(i) - ENDDO -! -! Iteration of the variables at the reference level zref : corrector calculation ; see Hess & McAvaney, 1995 -! - DO n = 1, niter -! - okri=.TRUE. - CALL screenc(klon, knon, nsrf, zxli, & - & u_zref, temp, q_zref, zref, & - & ts1, qsurf, z0m, z0h, psol, & - & ustar, testar, qstar, okri, ri1, & - & pref, delu, delte, delq) -! - DO i = 1, knon - u_zref(i) = delu(i) - q_zref(i) = delq(i) + max(qsurf(i),0.0) - te_zref(i) = delte(i) + ts1(i) -! -! return to normal temperature -! - temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA) -! temp(i) = te_zref(i) - (zref* RG)/RCPD/ & -! (1 + RVTMP2 * max(q_zref(i),0.0)) -! -!IM +++ -! IF(temp(i).GT.350.) THEN -! WRITE(*,*) 'temp(i) GT 350 K !!',i,nsrf,temp(i) -! ENDIF -!IM --- -! - IF(n.EQ.ncon) THEN - te_zref_con(i) = te_zref(i) - q_zref_con(i) = q_zref(i) - ENDIF -! - ENDDO -! - ENDDO -! -! verifier le critere de convergence : 0.25% pour te_zref et 5% pour qe_zref -! -! DO i = 1, knon -! conv_te(i) = (te_zref(i) - te_zref_con(i))/te_zref_con(i) -! conv_q(i) = (q_zref(i) - q_zref_con(i))/q_zref_con(i) -!IM +++ -! IF(abs(conv_te(i)).GE.0.0025.AND.abs(conv_q(i)).GE.0.05) THEN -! PRINT*,'DIV','i=',i,te_zref_con(i),te_zref(i),conv_te(i), & -! q_zref_con(i),q_zref(i),conv_q(i) -! ENDIF -!IM --- -! ENDDO -! - DO i = 1, knon - q_zref_c(i) = q_zref(i) - temp_c(i) = temp(i) -! -! IF(zri1(i).LT.0.) THEN -! IF(nsrf.EQ.1) THEN -! ok_pred(i)=1. -! ok_corr(i)=0. -! ELSE -! ok_pred(i)=0. -! ok_corr(i)=1. -! ENDIF -! ELSE -! ok_pred(i)=0. -! ok_corr(i)=1. -! ENDIF -! - ok_pred(i)=0. - ok_corr(i)=1. -! - t_2m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i) - q_2m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i) -!IM +++ -! IF(n.EQ.niter) THEN -! IF(t_2m(i).LT.t1(i).AND.t_2m(i).LT.ts1(i)) THEN -! PRINT*,' BAD t2m LT ',i,nsrf,t_2m(i),t1(i),ts1(i) -! ELSEIF(t_2m(i).GT.t1(i).AND.t_2m(i).GT.ts1(i)) THEN -! PRINT*,' BAD t2m GT ',i,nsrf,t_2m(i),t1(i),ts1(i) -! ENDIF -! ENDIF -!IM --- - ENDDO -! -! -!----------First aproximation of variables at zref -------------------------- -! - zref = 10.0 - CALL screenp(klon, knon, nsrf, speed, tpot, q1, & - & ts1, qsurf, z0m, lmon, & - & ustar, testar, qstar, zref, & - & delu, delte, delq) -! - DO i = 1, knon - u_zref(i) = delu(i) - q_zref(i) = max(qsurf(i),0.0) + delq(i) - te_zref(i) = ts1(i) + delte(i) - temp(i) = te_zref(i) * (psol(i)/pat1(i))**(-RKAPPA) -! temp(i) = te_zref(i) - (zref* RG)/RCPD/ & -! (1 + RVTMP2 * max(q_zref(i),0.0)) - u_zref_p(i) = u_zref(i) - ENDDO -! -! Iteration of the variables at the reference level zref : corrector ; see Hess & McAvaney, 1995 -! - DO n = 1, niter -! - okri=.TRUE. - CALL screenc(klon, knon, nsrf, zxli, & - & u_zref, temp, q_zref, zref, & - & ts1, qsurf, z0m, z0h, psol, & - & ustar, testar, qstar, okri, ri1, & - & pref, delu, delte, delq) -! - DO i = 1, knon - u_zref(i) = delu(i) - q_zref(i) = delq(i) + max(qsurf(i),0.0) - te_zref(i) = delte(i) + ts1(i) - temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA) -! temp(i) = te_zref(i) - (zref* RG)/RCPD/ & -! (1 + RVTMP2 * max(q_zref(i),0.0)) - ENDDO -! - ENDDO -! - DO i = 1, knon - u_zref_c(i) = u_zref(i) -! - u_10m(i) = u_zref_p(i) * ok_pred(i) + u_zref_c(i) * ok_corr(i) -! -!AM - q_zref_c(i) = q_zref(i) - temp_c(i) = temp(i) - t_10m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i) - q_10m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i) -!MA - ENDDO -! - RETURN - END subroutine stdlevvar