Index: LMDZ6/trunk/libf/phylmd/Dust/checknanqfi.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/Dust/checknanqfi.f90 (revision 5560) +++ LMDZ6/trunk/libf/phylmd/Dust/checknanqfi.f90 (revision 5561) @@ -1,5 +1,4 @@ SUBROUTINE checknanqfi(zq,qmin,qmax,comment) USE dimphy - USE, intrinsic :: ieee_arithmetic IMPLICIT NONE @@ -17,5 +16,5 @@ DO i = 1, klon ! IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN - IF (ieee_is_nan(zq(i,k))) THEN + IF (isnan(zq(i,k))) THEN jbad = jbad + 1 jadrs(jbad) = i Index: LMDZ6/trunk/libf/phylmd/acama_gwd_rando_m.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/acama_gwd_rando_m.f90 (revision 5560) +++ LMDZ6/trunk/libf/phylmd/acama_gwd_rando_m.f90 (revision 5561) @@ -5,23 +5,129 @@ USE clesphys_mod_h - IMPLICIT NONE - LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true. - LOGICAL, SAVE :: firstcall = .TRUE. + implicit none + +contains + + SUBROUTINE ACAMA_GWD_rando(DTIME, pp, plat, tt, uu, vv, rot, & + zustr, zvstr, d_u, d_v,east_gwstress,west_gwstress) + + ! Parametrization of the momentum flux deposition due to a discrete + ! number of gravity waves. + ! Author: F. Lott, A. de la Camara + ! July, 24th, 2014 + ! Gaussian distribution of the source, source is vorticity squared + ! Reference: de la Camara and Lott (GRL, 2015, vol 42, 2071-2078 ) + ! Lott et al (JAS, 2010, vol 67, page 157-170) + ! Lott et al (JAS, 2012, vol 69, page 2134-2151) + +! ONLINE: +USE yoegwd_mod_h + USE yomcst_mod_h +use dimphy, only: klon, klev + use assert_m, only: assert + USE ioipsl_getin_p_mod, ONLY : getin_p + USE vertical_layers_mod, ONLY : presnivs + + +! OFFLINE: +! include "dimensions_mod.f90" +! include "dimphy.h" +!END DIFFERENCE + + ! 0. DECLARATIONS: + + ! 0.1 INPUTS + REAL, intent(in)::DTIME ! Time step of the Physics + REAL, intent(in):: PP(:, :) ! (KLON, KLEV) Pressure at full levels + REAL, intent(in):: ROT(:,:) ! Relative vorticity + REAL, intent(in):: TT(:, :) ! (KLON, KLEV) Temp at full levels + REAL, intent(in):: UU(:, :) ! (KLON, KLEV) Zonal wind at full levels + REAL, intent(in):: VV(:, :) ! (KLON, KLEV) Merid wind at full levels + REAL, intent(in):: PLAT(:) ! (KLON) LATITUDE + + ! 0.2 OUTPUTS + REAL, intent(out):: zustr(:), zvstr(:) ! (KLON) Surface Stresses + + REAL, intent(inout):: d_u(:, :), d_v(:, :) + REAL, intent(inout):: east_gwstress(:, :) ! Profile of eastward stress + REAL, intent(inout):: west_gwstress(:, :) ! Profile of westward stress + ! (KLON, KLEV) tendencies on winds + + ! O.3 INTERNAL ARRAYS + REAL BVLOW(klon) ! LOW LEVEL BV FREQUENCY + REAL ROTBA(KLON),CORIO(KLON) ! BAROTROPIC REL. VORTICITY AND PLANETARY + REAL UZ(KLON, KLEV + 1) + + INTEGER II, JJ, LL + + ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED + + REAL DELTAT + + ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS + + INTEGER, PARAMETER:: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO + INTEGER JK, JP, JO, JW + INTEGER, PARAMETER:: NA = 5 !number of realizations to get the phase speed + REAL KMIN, KMAX ! Min and Max horizontal wavenumbers + REAL CMIN, CMAX ! Min and Max absolute ph. vel. + REAL CPHA ! absolute PHASE VELOCITY frequency + REAL ZK(NW, KLON) ! Horizontal wavenumber amplitude + REAL ZP(NW, KLON) ! Horizontal wavenumber angle + REAL ZO(NW, KLON) ! Absolute frequency ! + + ! Waves Intr. freq. at the 1/2 lev surrounding the full level + REAL ZOM(NW, KLON), ZOP(NW, KLON) + + ! Wave EP-fluxes at the 2 semi levels surrounding the full level + REAL WWM(NW, KLON), WWP(NW, KLON) + + REAL RUW0(NW, KLON) ! Fluxes at launching level + + REAL RUWP(NW, KLON), RVWP(NW, KLON) + ! Fluxes X and Y for each waves at 1/2 Levels + + INTEGER LAUNCH, LTROP ! Launching altitude and tropo altitude + + REAL XLAUNCH ! Controle the launching altitude + REAL XTROP ! SORT of Tropopause altitude + REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels + REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels + + REAL PRMAX ! Maximum value of PREC, and for which our linear formula + ! for GWs parameterisation apply + + ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS + + REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ + REAL CORSEC ! SECURITY FOR INTRINSIC CORIOLIS + REAL RUWFRT,SATFRT + + ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE + + REAL H0 ! Characteristic Height of the atmosphere + REAL DZ ! Characteristic depth of the source! + REAL PR, TR ! Reference Pressure and Temperature + + REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude + + REAL UH(KLON, KLEV + 1), VH(KLON, KLEV + 1) ! Winds at 1/2 levels + REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels + REAL PSEC ! Security to avoid division by 0 pressure + REAL PHM1(KLON, KLEV + 1) ! 1/Press at 1/2 levels + REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels + REAL BVSEC ! Security to avoid negative BVF + + REAL, DIMENSION(klev+1) ::HREF + LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true. + LOGICAL, SAVE :: firstcall = .TRUE. !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp) - - INTEGER, PARAMETER:: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO - INTEGER, PARAMETER:: NA = 5 !number of realizations to get the phase speed - - -contains - - SUBROUTINE ACAMA_GWD_rando_first - use dimphy, only: klev - USE ioipsl_getin_p_mod, ONLY : getin_p - IMPLICIT NONE - CHARACTER (LEN=20),PARAMETER :: modname='acama_gwd_rando_m' + + CHARACTER (LEN=20) :: modname='acama_gwd_rando_m' CHARACTER (LEN=80) :: abort_message - - IF (firstcall) THEN + + + + IF (firstcall) THEN ! Cle introduite pour resoudre un probleme de non reproductibilite ! Le but est de pouvoir tester de revenir a la version precedenete @@ -34,119 +140,5 @@ firstcall=.false. ! CALL iophys_ini(dtime) - ENDIF - END SUBROUTINE ACAMA_GWD_rando_first - - SUBROUTINE ACAMA_GWD_rando(DTIME, pp, plat, tt, uu, vv, rot, & - zustr, zvstr, d_u, d_v,east_gwstress,west_gwstress) - - ! Parametrization of the momentum flux deposition due to a discrete - ! number of gravity waves. - ! Author: F. Lott, A. de la Camara - ! July, 24th, 2014 - ! Gaussian distribution of the source, source is vorticity squared - ! Reference: de la Camara and Lott (GRL, 2015, vol 42, 2071-2078 ) - ! Lott et al (JAS, 2010, vol 67, page 157-170) - ! Lott et al (JAS, 2012, vol 69, page 2134-2151) - -! ONLINE: -USE yoegwd_mod_h - USE yomcst_mod_h -use dimphy, only: klon, klev - use assert_m, only: assert - USE vertical_layers_mod, ONLY : presnivs - USE lmdz_fake_call, ONLY : fake_call - -! OFFLINE: -! include "dimensions_mod.f90" -! include "dimphy.h" -!END DIFFERENCE - - ! 0. DECLARATIONS: - - ! 0.1 INPUTS - REAL, intent(in)::DTIME ! Time step of the Physics - REAL, intent(in):: PP(KLON, KLEV) ! (KLON, KLEV) Pressure at full levels - REAL, intent(in):: ROT(KLON,KLEV) ! Relative vorticity - REAL, intent(in):: TT(KLON, KLEV) ! (KLON, KLEV) Temp at full levels - REAL, intent(in):: UU(KLON, KLEV) ! (KLON, KLEV) Zonal wind at full levels - REAL, intent(in):: VV(KLON, KLEV) ! (KLON, KLEV) Merid wind at full levels - REAL, intent(in):: PLAT(KLON) ! (KLON) LATITUDE - - ! 0.2 OUTPUTS - REAL, intent(out):: zustr(KLON), zvstr(KLON) ! (KLON) Surface Stresses - - REAL, intent(inout):: d_u(KLON, KLEV), d_v(KLON, KLEV) - REAL, intent(inout):: east_gwstress(KLON, KLEV) ! Profile of eastward stress - REAL, intent(inout):: west_gwstress(KLON, KLEV) ! Profile of westward stress - ! (KLON, KLEV) tendencies on winds - - ! O.3 INTERNAL ARRAYS - REAL BVLOW(klon) ! LOW LEVEL BV FREQUENCY - REAL ROTBA(KLON),CORIO(KLON) ! BAROTROPIC REL. VORTICITY AND PLANETARY - REAL UZ(KLON, KLEV + 1) - - INTEGER II, JJ, LL - - ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED - - REAL DELTAT - - ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS - - INTEGER JK, JP, JO, JW - REAL KMIN, KMAX ! Min and Max horizontal wavenumbers - REAL CMIN, CMAX ! Min and Max absolute ph. vel. - REAL CPHA ! absolute PHASE VELOCITY frequency - REAL ZK(KLON, NW) ! Horizontal wavenumber amplitude - REAL ZP(KLON, NW) ! Horizontal wavenumber angle - REAL ZO(KLON,NW) ! Absolute frequency ! - - ! Waves Intr. freq. at the 1/2 lev surrounding the full level - REAL ZOM(KLON, NW), ZOP(KLON, NW) - - ! Wave EP-fluxes at the 2 semi levels surrounding the full level - REAL WWM(KLON, NW), WWP(KLON, NW) - - REAL RUW0(KLON, NW) ! Fluxes at launching level - - REAL RUWP(KLON, NW), RVWP(KLON, NW) - ! Fluxes X and Y for each waves at 1/2 Levels - - INTEGER LAUNCH, LTROP ! Launching altitude and tropo altitude - - REAL XLAUNCH ! Controle the launching altitude - REAL XTROP ! SORT of Tropopause altitude - REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels - REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels - - REAL PRMAX ! Maximum value of PREC, and for which our linear formula - ! for GWs parameterisation apply - - ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS - - REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ - REAL CORSEC ! SECURITY FOR INTRINSIC CORIOLIS - REAL RUWFRT,SATFRT - - ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE - - REAL H0 ! Characteristic Height of the atmosphere - REAL DZ ! Characteristic depth of the source! - REAL PR, TR ! Reference Pressure and Temperature - - REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude - - REAL UH(KLON, KLEV + 1), VH(KLON, KLEV + 1) ! Winds at 1/2 levels - REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels - REAL PSEC ! Security to avoid division by 0 pressure - REAL PHM1(KLON, KLEV + 1) ! 1/Press at 1/2 levels - REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels - REAL BVSEC ! Security to avoid negative BVF - - REAL, DIMENSION(klev+1) ::HREF - - CHARACTER (LEN=20),PARAMETER :: modname='acama_gwd_rando_m' - CHARACTER (LEN=80) :: abort_message - + ENDIF !----------------------------------------------------------------- @@ -219,8 +211,4 @@ ! END ONLINE - CALL FAKE_CALL(BVLOW) ! to be suppress in future - CALL FAKE_CALL(CORIO) ! to be suppress in future - CALL FAKE_CALL(ROTBA) ! to be suppress in future - IF(DELTAT < DTIME)THEN ! PRINT *, 'flott_gwd_rando: deltat < dtime!' @@ -288,5 +276,5 @@ - TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1)) * RD / H0 end DO - BVLOW(:) = 0.5 * (TT(:, LTROP )+ TT(:, LAUNCH)) & + BVLOW = 0.5 * (TT(:, LTROP )+ TT(:, LAUNCH)) & * RD**2 / RCPD / H0**2 + (TT(:, LTROP ) & - TT(:, LAUNCH))/(ZH(:, LTROP )- ZH(:, LAUNCH)) * RD / H0 @@ -357,12 +345,12 @@ DO II = 1, KLON ! Angle (0 or PI so far) - ! ZP(II,JW) = (SIGN(1., 0.5 - MOD(TT(II, JW) * 10., 1.)) + 1.) & + ! ZP(JW, II) = (SIGN(1., 0.5 - MOD(TT(II, JW) * 10., 1.)) + 1.) & ! * RPI / 2. ! Angle between 0 and pi - ZP(II,JW) = MOD(TT(II, JW) * 10., 1.) * RPI + ZP(JW, II) = MOD(TT(II, JW) * 10., 1.) * RPI ! TEST WITH POSITIVE WAVES ONLY (Part I/II) -! ZP(II,JW) = 0. +! ZP(JW, II) = 0. ! Horizontal wavenumber amplitude - ZK(II,JW) = KMIN + (KMAX - KMIN) * MOD(TT(II, JW) * 100., 1.) + ZK(JW, II) = KMIN + (KMAX - KMIN) * MOD(TT(II, JW) * 100., 1.) ! Horizontal phase speed CPHA = 0. @@ -373,27 +361,27 @@ IF (CPHA.LT.0.) THEN CPHA = -1.*CPHA - ZP(II,JW) = ZP(II,JW) + RPI + ZP(JW,II) = ZP(JW,II) + RPI ! TEST WITH POSITIVE WAVES ONLY (Part II/II) -! ZP(II,JW) = 0. +! ZP(JW, II) = 0. ENDIF CPHA = CPHA + CMIN !we dont allow |c|<1m/s ! Absolute frequency is imposed - ZO(II, JW) = CPHA * ZK(II,JW) + ZO(JW, II) = CPHA * ZK(JW, II) ! Intrinsic frequency is imposed - ZO(II, JW) = ZO(II, JW) & - + ZK(II,JW) * COS(ZP(II,JW)) * UH(II, LAUNCH) & - + ZK(II,JW) * SIN(ZP(II,JW)) * VH(II, LAUNCH) + ZO(JW, II) = ZO(JW, II) & + + ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH) & + + ZK(JW, II) * SIN(ZP(JW, II)) * VH(II, LAUNCH) ! Momentum flux at launch lev ! LAUNCHED RANDOM WAVES WITH LOG-NORMAL AMPLITUDE ! RIGHT IN THE SH (GWD4 after 1990) - RUW0(II, JW) = 0. + RUW0(JW, II) = 0. DO JJ = 1, NA - RUW0(II, JW) = RUW0(II, JW) + & + RUW0(JW, II) = RUW0(JW,II) + & 2.*(MOD(TT(II, JW+4*(JJ-1)+JJ)**2, 1.)-0.5)*SQRT(3.)/SQRT(NA*1.) END DO - RUW0(II,JW) = RUWFRT & - * EXP(RUW0(II,JW))/1250. & ! 2 mpa at south pole + RUW0(JW, II) = RUWFRT & + * EXP(RUW0(JW,II))/1250. & ! 2 mpa at south pole *((1.05+SIN(PLAT(II)*RPI/180.))/(1.01+SIN(PLAT(II)*RPI/180.))-2.05/2.01) - ! RUW0(II,JW) = RUWFRT + ! RUW0(JW, II) = RUWFRT ENDDO end DO @@ -409,7 +397,7 @@ ! Evaluate intrinsic frequency at launching altitude: - ZOP(:, JW) = ZO(:, JW) & - - ZK(:,JW) * COS(ZP(:,JW)) * UH(:, LAUNCH) & - - ZK(:,JW) * SIN(ZP(:,JW)) * VH(:, LAUNCH) + ZOP(JW, :) = ZO(JW, :) & + - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) & + - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH) ! VERSION WITH FRONTAL SOURCES @@ -418,6 +406,6 @@ ! tanh limitation for values above CORIO (inertial instability). - ! WWP(:,JW) = RUW0(:,JW) & - WWP(:,JW) = RUWFRT & + ! WWP(JW, :) = RUW0(JW, :) & + WWP(JW, :) = RUWFRT & ! * (CORIO(:)*TANH(ROTBA(:)/CORIO(:)))**2 & ! * ABS((CORIO(:)*TANH(ROTBA(:)/CORIO(:)))*CORIO(:)) & @@ -425,6 +413,6 @@ ! * (CORIO(:)*CORIO(:)) & ! MODERATION BY THE DEPTH OF THE SOURCE (DZ HERE) - ! *EXP(-BVLOW(:)**2/MAX(ABS(ZOP(:,JW)),ZOISEC)**2 & - ! *ZK(:,JW)**2*DZ**2) & + ! *EXP(-BVLOW(:)**2/MAX(ABS(ZOP(JW, :)),ZOISEC)**2 & + ! *ZK(JW, :)**2*DZ**2) & ! COMPLETE FORMULA: !* CORIO(:)**2*TANH(ROTBA(:)/CORIO(:)**2) & @@ -432,5 +420,5 @@ ! RESTORE DIMENSION OF A FLUX ! *RD*TR/PR - ! *1. + RUW0(:,JW) + ! *1. + RUW0(JW, :) *1. @@ -441,6 +429,6 @@ ! Put the stress in the right direction: - RUWP(:,JW) = SIGN(1., ZOP(:,JW))*COS(ZP(:,JW)) * WWP(:,JW) - RVWP(:,JW) = SIGN(1., ZOP(:,JW))*SIN(ZP(:,JW)) * WWP(:,JW) + RUWP(JW, :) = SIGN(1., ZOP(JW, :))*COS(ZP(JW, :)) * WWP(JW, :) + RVWP(JW, :) = SIGN(1., ZOP(JW, :))*SIN(ZP(JW, :)) * WWP(JW, :) end DO @@ -452,6 +440,6 @@ RVW(:, LL) = 0 DO JW = 1, NW - RUW(:, LL) = RUW(:, LL) + RUWP(:,JW) - RVW(:, LL) = RVW(:, LL) + RVWP(:,JW) + RUW(:, LL) = RUW(:, LL) + RUWP(JW, :) + RVW(:, LL) = RVW(:, LL) + RVWP(JW, :) end DO end DO @@ -465,26 +453,26 @@ ! to the next) DO JW = 1, NW - ZOM(:, JW) = ZOP(:,JW) - WWM(:,JW) = WWP(:,JW) + ZOM(JW, :) = ZOP(JW, :) + WWM(JW, :) = WWP(JW, :) ! Intrinsic Frequency - ZOP(:,JW) = ZO(:, JW) - ZK(:,JW) * COS(ZP(:,JW)) * UH(:, LL + 1) & - - ZK(:,JW) * SIN(ZP(:,JW)) * VH(:, LL + 1) + ZOP(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LL + 1) & + - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LL + 1) ! No breaking (Eq.6) ! Dissipation (Eq. 8) - WWP(:,JW) = WWM(:,JW) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) & + WWP(JW, :) = WWM(JW, :) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) & + PH(:, LL)) * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 & - / MAX(ABS(ZOP(:,JW) + ZOM(:, JW)) / 2., ZOISEC)**4 & - * ZK(:,JW)**3 * (ZH(:, LL + 1) - ZH(:, LL))) + / MAX(ABS(ZOP(JW, :) + ZOM(JW, :)) / 2., ZOISEC)**4 & + * ZK(JW, :)**3 * (ZH(:, LL + 1) - ZH(:, LL))) ! Critical levels (forced to zero if intrinsic frequency changes sign) ! Saturation (Eq. 12) - WWP(:,JW) = min(WWP(:,JW), MAX(0., & - SIGN(1., ZOP(:,JW) * ZOM(:, JW))) * ABS(ZOP(:,JW))**3 & + WWP(JW, :) = min(WWP(JW, :), MAX(0., & + SIGN(1., ZOP(JW, :) * ZOM(JW, :))) * ABS(ZOP(JW, :))**3 & ! / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * SATFRT**2 * KMIN**2 & / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * KMIN**2 & ! *(SATFRT*(2.5+1.5*TANH((ZH(:,LL+1)/H0-8.)/2.)))**2 & *SATFRT**2 & - / ZK(:,JW)**4) + / ZK(JW, :)**4) end DO @@ -493,6 +481,6 @@ DO JW = 1, NW - RUWP(:,JW) = SIGN(1., ZOP(:,JW))*COS(ZP(:,JW)) * WWP(:,JW) - RVWP(:,JW) = SIGN(1., ZOP(:,JW))*SIN(ZP(:,JW)) * WWP(:,JW) + RUWP(JW, :) = SIGN(1., ZOP(JW, :))*COS(ZP(JW, :)) * WWP(JW, :) + RVWP(JW, :) = SIGN(1., ZOP(JW, :))*SIN(ZP(JW, :)) * WWP(JW, :) end DO @@ -501,8 +489,8 @@ DO JW = 1, NW - RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(:,JW) - RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(:,JW) - EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,RUWP(:,JW))/REAL(NW) - WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,RUWP(:,JW))/REAL(NW) + RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(JW, :) + RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(JW, :) + EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,RUWP(JW,:))/FLOAT(NW) + WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,RUWP(JW,:))/FLOAT(NW) end DO end DO Index: LMDZ6/trunk/libf/phylmd/flott_gwd_rando_m.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/flott_gwd_rando_m.f90 (revision 5560) +++ LMDZ6/trunk/libf/phylmd/flott_gwd_rando_m.f90 (revision 5561) @@ -6,20 +6,121 @@ USE clesphys_mod_h implicit none + +contains + + SUBROUTINE FLOTT_GWD_rando(DTIME, pp, tt, uu, vv, prec, zustr, zvstr, d_u, & + d_v,east_gwstress,west_gwstress) + + ! Parametrization of the momentum flux deposition due to a discrete + ! number of gravity waves. + ! Author: F. Lott + ! July, 12th, 2012 + ! Gaussian distribution of the source, source is precipitation + ! Reference: Lott (JGR, vol 118, page 8897, 2013) + + !ONLINE: + USE yomcst_mod_h +use dimphy, only: klon, klev + use assert_m, only: assert + USE ioipsl_getin_p_mod, ONLY : getin_p + USE vertical_layers_mod, ONLY : presnivs + USE yoegwd_mod_h + CHARACTER (LEN=20) :: modname='flott_gwd_rando' + CHARACTER (LEN=80) :: abort_message + + ! OFFLINE: + ! include "dimensions_mod.f90" + ! include "dimphy.h" + ! END OF DIFFERENCE ONLINE-OFFLINE + + ! 0. DECLARATIONS: + + ! 0.1 INPUTS + REAL, intent(in)::DTIME ! Time step of the Physics + REAL, intent(in):: pp(:, :) ! (KLON, KLEV) Pressure at full levels + REAL, intent(in):: prec(:) ! (klon) Precipitation (kg/m^2/s) + REAL, intent(in):: TT(:, :) ! (KLON, KLEV) Temp at full levels + REAL, intent(in):: UU(:, :) ! (KLON, KLEV) Zonal wind at full levels + REAL, intent(in):: VV(:, :) ! (KLON, KLEV) Merid wind at full levels + + ! 0.2 OUTPUTS + REAL, intent(out):: zustr(:), zvstr(:) ! (KLON) Surface Stresses + + REAL, intent(inout):: d_u(:, :), d_v(:, :) + REAL, intent(inout):: east_gwstress(:, :) ! Profile of eastward stress + REAL, intent(inout):: west_gwstress(:, :) ! Profile of westward stress + + ! (KLON, KLEV) tendencies on winds + + ! O.3 INTERNAL ARRAYS + REAL BVLOW(klon) + REAL DZ ! Characteristic depth of the Source + + INTEGER II, JJ, LL + + ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED + + REAL DELTAT + + ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS + INTEGER, PARAMETER:: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO + INTEGER JK, JP, JO, JW INTEGER, PARAMETER:: NA = 5 !number of realizations to get the phase speed + REAL KMIN, KMAX ! Min and Max horizontal wavenumbers + REAL CMAX ! standard deviation of the phase speed distribution + REAL RUWMAX,SAT ! ONLINE SPECIFIED IN run.def + REAL CPHA ! absolute PHASE VELOCITY frequency + REAL ZK(NW, KLON) ! Horizontal wavenumber amplitude + REAL ZP(NW, KLON) ! Horizontal wavenumber angle + REAL ZO(NW, KLON) ! Absolute frequency ! + + ! Waves Intr. freq. at the 1/2 lev surrounding the full level + REAL ZOM(NW, KLON), ZOP(NW, KLON) + + ! Wave EP-fluxes at the 2 semi levels surrounding the full level + REAL WWM(NW, KLON), WWP(NW, KLON) + + REAL RUW0(NW, KLON) ! Fluxes at launching level + + REAL RUWP(NW, KLON), RVWP(NW, KLON) + ! Fluxes X and Y for each waves at 1/2 Levels + + INTEGER LAUNCH, LTROP ! Launching altitude and tropo altitude + + REAL XLAUNCH ! Controle the launching altitude + REAL XTROP ! SORT of Tropopause altitude + REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels + REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels + + REAL PRMAX ! Maximum value of PREC, and for which our linear formula + ! for GWs parameterisation apply + + ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS + + REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ + + ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE + + REAL H0 ! Characteristic Height of the atmosphere + REAL PR, TR ! Reference Pressure and Temperature + + REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude + + REAL UH(KLON, KLEV + 1), VH(KLON, KLEV + 1) ! Winds at 1/2 levels + REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels + REAL PSEC ! Security to avoid division by 0 pressure + REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels + REAL BVSEC ! Security to avoid negative BVF + REAL RAN_NUM_1,RAN_NUM_2,RAN_NUM_3 + + REAL, DIMENSION(klev+1) ::HREF + LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true. LOGICAL, SAVE :: firstcall = .TRUE. - !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp) - -contains - - SUBROUTINE FLOTT_GWD_rando_first - use dimphy, only: klev - USE ioipsl_getin_p_mod, ONLY : getin_p - IMPLICIT NONE - CHARACTER (LEN=20),PARAMETER :: modname='acama_gwd_rando_m' - CHARACTER (LEN=80) :: abort_message - - IF (firstcall) THEN + !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp) + + + IF (firstcall) THEN ! Cle introduite pour resoudre un probleme de non reproductibilite ! Le but est de pouvoir tester de revenir a la version precedenete @@ -32,114 +133,4 @@ firstcall=.false. ENDIF - END SUBROUTINE FLOTT_GWD_rando_first - - - SUBROUTINE FLOTT_GWD_rando(DTIME, PP, tt, uu, vv, prec, zustr, zvstr, d_u, & - d_v,east_gwstress,west_gwstress) - - ! Parametrization of the momentum flux deposition due to a discrete - ! number of gravity waves. - ! Author: F. Lott - ! July, 12th, 2012 - ! Gaussian distribution of the source, source is precipitation - ! Reference: Lott (JGR, vol 118, page 8897, 2013) - - !ONLINE: - USE yomcst_mod_h -use dimphy, only: klon, klev - use assert_m, only: assert - USE ioipsl_getin_p_mod, ONLY : getin_p - USE vertical_layers_mod, ONLY : presnivs - USE yoegwd_mod_h - USE lmdz_fake_call, ONLY : fake_call - - CHARACTER (LEN=20),PARAMETER :: modname='flott_gwd_rando' - CHARACTER (LEN=80) :: abort_message - - ! OFFLINE: - ! include "dimensions_mod.f90" - ! include "dimphy.h" - ! END OF DIFFERENCE ONLINE-OFFLINE - - ! 0. DECLARATIONS: - - ! 0.1 INPUTS - REAL, intent(in)::DTIME ! Time step of the Physics - REAL, intent(in):: pp(KLON, KLEV) ! (KLON, KLEV) Pressure at full levels - REAL, intent(in):: prec(KLON) ! (klon) Precipitation (kg/m^2/s) - REAL, intent(in):: TT(KLON, KLEV) ! (KLON, KLEV) Temp at full levels - REAL, intent(in):: UU(KLON, KLEV) ! (KLON, KLEV) Zonal wind at full levels - REAL, intent(in):: VV(KLON, KLEV) ! (KLON, KLEV) Merid wind at full levels - - ! 0.2 OUTPUTS - REAL, intent(out):: zustr(KLON), zvstr(KLON) ! (KLON) Surface Stresses - - REAL, intent(inout):: d_u(KLON, KLEV), d_v(KLON, KLEV) - REAL, intent(inout):: east_gwstress(KLON, KLEV) ! Profile of eastward stress - REAL, intent(inout):: west_gwstress(KLON, KLEV) ! Profile of westward stress - - ! (KLON, KLEV) tendencies on winds - - ! O.3 INTERNAL ARRAYS - REAL BVLOW(klon) - REAL DZ ! Characteristic depth of the Source - - INTEGER II, JJ, LL - - ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED - - REAL DELTAT - - ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS - - INTEGER JK, JP, JO, JW - REAL KMIN, KMAX ! Min and Max horizontal wavenumbers - REAL CMAX ! standard deviation of the phase speed distribution - REAL RUWMAX,SAT ! ONLINE SPECIFIED IN run.def - REAL CPHA ! absolute PHASE VELOCITY frequency - REAL ZK(KLON, NW) ! Horizontal wavenumber amplitude - REAL ZP(KLON, NW) ! Horizontal wavenumber angle - REAL ZO(KLON, NW) ! Absolute frequency ! - - ! Waves Intr. freq. at the 1/2 lev surrounding the full level - REAL ZOM(KLON, NW), ZOP(KLON, NW) - - ! Wave EP-fluxes at the 2 semi levels surrounding the full level - REAL WWM(KLON, NW), WWP(KLON, NW) - - REAL RUW0(KLON, NW) ! Fluxes at launching level - - REAL RUWP(KLON, NW), RVWP(KLON, NW) - ! Fluxes X and Y for each waves at 1/2 Levels - - INTEGER LAUNCH, LTROP ! Launching altitude and tropo altitude - - REAL XLAUNCH ! Controle the launching altitude - REAL XTROP ! SORT of Tropopause altitude - REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels - REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels - - REAL PRMAX ! Maximum value of PREC, and for which our linear formula - ! for GWs parameterisation apply - - ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS - - REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ - - ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE - - REAL H0 ! Characteristic Height of the atmosphere - REAL PR, TR ! Reference Pressure and Temperature - - REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude - - REAL UH(KLON, KLEV + 1), VH(KLON, KLEV + 1) ! Winds at 1/2 levels - REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels - REAL PSEC ! Security to avoid division by 0 pressure - REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels - REAL BVSEC ! Security to avoid negative BVF - REAL RAN_NUM_1,RAN_NUM_2,RAN_NUM_3 - - REAL, DIMENSION(klev+1) ::HREF @@ -206,5 +197,5 @@ !END ONLINE ENDIF - + IF(DELTAT < DTIME)THEN abort_message='flott_gwd_rando: deltat < dtime!' @@ -216,7 +207,5 @@ CALL abort_physic(modname,abort_message,1) ENDIF - - CALL FAKE_CALL(BVLOW) ! to be suppress in future - + ! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS @@ -303,8 +292,8 @@ RAN_NUM_1=MOD(TT(II, JW) * 10., 1.) RAN_NUM_2= MOD(TT(II, JW) * 100., 1.) - ZP(II, JW) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.) & + ZP(JW, II) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.) & * RPI / 2. ! Horizontal wavenumber amplitude - ZK(II, JW) = KMIN + (KMAX - KMIN) *RAN_NUM_2 + ZK(JW, II) = KMIN + (KMAX - KMIN) *RAN_NUM_2 ! Horizontal phase speed CPHA = 0. @@ -316,14 +305,14 @@ IF (CPHA.LT.0.) THEN CPHA = -1.*CPHA - ZP(II, JW) = ZP(II, JW) + RPI + ZP(JW,II) = ZP(JW,II) + RPI ENDIF ! Absolute frequency is imposed - ZO(II, JW) = CPHA * ZK(II, JW) + ZO(JW, II) = CPHA * ZK(JW, II) ! Intrinsic frequency is imposed - ZO(II, JW) = ZO(II, JW) & - + ZK(II, JW) * COS(ZP(II, JW)) * UH(II, LAUNCH) & - + ZK(II, JW) * SIN(ZP(II, JW)) * VH(II, LAUNCH) + ZO(JW, II) = ZO(JW, II) & + + ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH) & + + ZK(JW, II) * SIN(ZP(JW, II)) * VH(II, LAUNCH) ! Momentum flux at launch lev - RUW0(II, JW) = RUWMAX + RUW0(JW, II) = RUWMAX ENDDO ENDDO @@ -337,7 +326,7 @@ ! Evaluate intrinsic frequency at launching altitude: - ZOP(:,JW) = ZO(:, JW) & - - ZK(:, JW) * COS(ZP(:, JW)) * UH(:, LAUNCH) & - - ZK(:, JW) * SIN(ZP(:, JW)) * VH(:, LAUNCH) + ZOP(JW, :) = ZO(JW, :) & + - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) & + - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH) ! VERSION WITH CONVECTIVE SOURCE @@ -348,21 +337,21 @@ ! tanh limitation to values above prmax: - WWP(:, JW) = RUW0(:, JW) & + WWP(JW, :) = RUW0(JW, :) & * (RD / RCPD / H0 * RLVTT * PRMAX * TANH(PREC(:) / PRMAX))**2 ! Factor related to the characteristics of the waves: - WWP(:, JW) = WWP(:, JW) * ZK(:, JW)**3 / KMIN / BVLOW(:) & - / MAX(ABS(ZOP(:, JW)), ZOISEC)**3 + WWP(JW, :) = WWP(JW, :) * ZK(JW, :)**3 / KMIN / BVLOW(:) & + / MAX(ABS(ZOP(JW, :)), ZOISEC)**3 ! Moderation by the depth of the source (dz here): - WWP(:, JW) = WWP(:, JW) & - * EXP(- BVLOW(:)**2 / MAX(ABS(ZOP(:, JW)), ZOISEC)**2 * ZK(:, JW)**2 & + WWP(JW, :) = WWP(JW, :) & + * EXP(- BVLOW(:)**2 / MAX(ABS(ZOP(JW, :)), ZOISEC)**2 * ZK(JW, :)**2 & * DZ**2) ! Put the stress in the right direction: - RUWP(:, JW) = ZOP(:, JW) / MAX(ABS(ZOP(:, JW)), ZOISEC)**2 & - * BV(:, LAUNCH) * COS(ZP(:, JW)) * WWP(:, JW)**2 - RVWP(:, JW) = ZOP(:, JW) / MAX(ABS(ZOP(:, JW)), ZOISEC)**2 & - * BV(:, LAUNCH) * SIN(ZP(:, JW)) * WWP(:, JW)**2 + RUWP(JW, :) = ZOP(JW, :) / MAX(ABS(ZOP(JW, :)), ZOISEC)**2 & + * BV(:, LAUNCH) * COS(ZP(JW, :)) * WWP(JW, :)**2 + RVWP(JW, :) = ZOP(JW, :) / MAX(ABS(ZOP(JW, :)), ZOISEC)**2 & + * BV(:, LAUNCH) * SIN(ZP(JW, :)) * WWP(JW, :)**2 end DO @@ -374,6 +363,6 @@ RVW(:, LL) = 0 DO JW = 1, NW - RUW(:, LL) = RUW(:, LL) + RUWP(:, JW) - RVW(:, LL) = RVW(:, LL) + RVWP(:, JW) + RUW(:, LL) = RUW(:, LL) + RUWP(JW, :) + RVW(:, LL) = RVW(:, LL) + RVWP(JW, :) end DO end DO @@ -387,23 +376,23 @@ ! to the next) DO JW = 1, NW - ZOM(:, JW) = ZOP(:,JW) - WWM(:, JW) = WWP(:, JW) + ZOM(JW, :) = ZOP(JW, :) + WWM(JW, :) = WWP(JW, :) ! Intrinsic Frequency - ZOP(:, JW) = ZO(:, JW) - ZK(:, JW) * COS(ZP(:, JW)) * UH(:, LL + 1) & - - ZK(:, JW) * SIN(ZP(:, JW)) * VH(:, LL + 1) + ZOP(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LL + 1) & + - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LL + 1) ! No breaking (Eq.6) ! Dissipation (Eq. 8) - WWP(:, JW) = WWM(:, JW) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) & + WWP(JW, :) = WWM(JW, :) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) & + PH(:, LL)) * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 & - / MAX(ABS(ZOP(:, JW) + ZOM(:, JW)) / 2., ZOISEC)**4 & - * ZK(:, JW)**3 * (ZH(:, LL + 1) - ZH(:, LL))) + / MAX(ABS(ZOP(JW, :) + ZOM(JW, :)) / 2., ZOISEC)**4 & + * ZK(JW, :)**3 * (ZH(:, LL + 1) - ZH(:, LL))) ! Critical levels (forced to zero if intrinsic frequency changes sign) ! Saturation (Eq. 12) - WWP(:, JW) = min(WWP(:, JW), MAX(0., & - SIGN(1., ZOP(:, JW) * ZOM(:, JW))) * ABS(ZOP(:, JW))**3 & + WWP(JW, :) = min(WWP(JW, :), MAX(0., & + SIGN(1., ZOP(JW, :) * ZOM(JW, :))) * ABS(ZOP(JW, :))**3 & / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * KMIN**2 & - * SAT**2 / ZK(:, JW)**4) + * SAT**2 / ZK(JW, :)**4) end DO @@ -412,6 +401,6 @@ DO JW = 1, NW - RUWP(:, JW) = SIGN(1., ZOP(:, JW))*COS(ZP(:, JW)) * WWP(:, JW) - RVWP(:, JW) = SIGN(1., ZOP(:, JW))*SIN(ZP(:, JW)) * WWP(:, JW) + RUWP(JW, :) = SIGN(1., ZOP(JW, :))*COS(ZP(JW, :)) * WWP(JW, :) + RVWP(JW, :) = SIGN(1., ZOP(JW, :))*SIN(ZP(JW, :)) * WWP(JW, :) end DO @@ -420,8 +409,8 @@ DO JW = 1, NW - RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(:, JW) - RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(:, JW) - EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,RUWP(:, JW))/REAL(NW) - WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,RUWP(:, JW))/REAL(NW) + RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(JW, :) + RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(JW, :) + EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,RUWP(JW,:))/FLOAT(NW) + WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,RUWP(JW,:))/FLOAT(NW) end DO end DO Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 5560) +++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 5561) @@ -20,5 +20,5 @@ ! PLEASE try to follow this rule - USE ACAMA_GWD_rando_m, only: ACAMA_GWD_rando, ACAMA_GWD_rando_first + USE ACAMA_GWD_rando_m, only: ACAMA_GWD_rando USE aero_mod USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, & @@ -32,5 +32,5 @@ USE dimphy USE etat0_limit_unstruct_mod - USE FLOTT_GWD_rando_m, only: FLOTT_GWD_rando, FLOTT_GWD_rando_first + USE FLOTT_GWD_rando_m, only: FLOTT_GWD_rando USE fonte_neige_mod, ONLY : fonte_neige_get_vars USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg @@ -4970,5 +4970,4 @@ IF (.not. ok_hines .and. ok_gwd_rando) then ! ym missing init for east_gwstress & west_gwstress -> added in phys_local_var_mod - CALL acama_GWD_rando_first() CALL acama_GWD_rando(PHYS_TSTEP, pplay, latitude_deg, t_seri, u_seri, & v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, & @@ -4989,5 +4988,4 @@ IF (ok_gwd_rando) THEN - CALL FLOTT_GWD_rando_first() CALL FLOTT_GWD_rando(PHYS_TSTEP, pplay, t_seri, u_seri, v_seri, & rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &