Index: /trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90 (revision 2298) +++ /trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90 (revision 2299) @@ -61,4 +61,7 @@ logical,save :: haze !$OMP THREADPRIVATE(photochem) + logical,save :: calllott_nonoro + logical,save :: gwd_convective_source +!$OMP THREADPRIVATE(calllott_nonoro,gwd_convective_source) integer,save :: iddist @@ -135,4 +138,6 @@ real,save :: surfemis !$OMP THREADPRIVATE(surfalbedo,surfemis) + real,save :: epflux_max +!$OMP THREADPRIVATE(epflux_max) logical,save :: iscallphys=.false.!existence of callphys.def Index: /trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90 (revision 2298) +++ /trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90 (revision 2299) @@ -22,4 +22,5 @@ use planetwide_mod, only: planetwide_sumval use callkeys_mod + use ioipsl_getin_p_mod, only : getin_p use mod_phys_lmdz_para, only : is_parallel @@ -290,4 +291,14 @@ write(*,*) " calladj = ",calladj + write(*,*)"call Lott's non-oro GWs parameterisation scheme ?" + calllott_nonoro=.false. ! default value + call getin_p("calllott_nonoro",calllott_nonoro) + write(*,*)" calllott_nonoro = ",calllott_nonoro + + write(*,*)"Max Value of EP flux at launching altitude" + epflux_max=0.0 ! default value + call getin_p("epflux_max",epflux_max) + write(*,*)"epflux_max = ",epflux_max + write(*,*) "call thermal plume model ?" calltherm=.false. ! default value Index: /trunk/LMDZ.GENERIC/libf/phystd/nonoro_gwd_ran_mod.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/nonoro_gwd_ran_mod.F90 (revision 2299) +++ /trunk/LMDZ.GENERIC/libf/phystd/nonoro_gwd_ran_mod.F90 (revision 2299) @@ -0,0 +1,491 @@ +MODULE nonoro_gwd_ran_mod + +IMPLICIT NONE + + REAL, allocatable, save :: du_nonoro_gwd(:, :) ! Zonal wind tendency due to GWD + REAL, allocatable, save :: dv_nonoro_gwd(:, :) ! Meridional wind tendency due to GWD + REAL, allocatable, save :: east_gwstress(:, :) ! Eastward stress profile + REAL, allocatable, save :: west_gwstress(:, :) ! Westward stress profile + +CONTAINS + + SUBROUTINE nonoro_gwd_ran(ngrid, nlayer, dtime, pp, & + zmax_therm, pt, pu, pv, pdt, pdu, pdv, & + zustr, zvstr, d_t, d_u, d_v) + +!-------------------------------------------------------------------------------- +! Parametrization of the momentum flux deposition due to a discrete +! number of gravity waves. +! F. Lott +! Version 14, Gaussian distribution of the source +! LMDz model online version +! ADAPTED FOR VENUS / F. LOTT + S. LEBONNOIS +! Version adapted on 03/04/2013: +! - input flux compensated in the deepest layers +! +! ADAPTED FOR MARS G.GILLI 02/2016 +! Revision with F.Forget 06/2016 Variable EP-flux +! according to +! PBL variation (max +! velocity thermals) +! D.BARDET 01/2020 - reproductibility of +! the launching altitude +! calculation +! - wave characteristic +! calculation using MOD +! - adding east_gwstress +! and west_gwstress variables +! +! ADAPTED FOR GENERIC D.BARDET 01/2020 +!--------------------------------------------------------------------------------- + + + + use comcstfi_mod, only: g, pi, cpp, r + USE ioipsl_getin_p_mod, ONLY : getin_p + use assert_m, only : assert + use vertical_layers_mod, only : presnivs + use callkeys_mod, only : epflux_max, & ! Max EP flux value at launching altitude (previous name: RUWMAX) + gwd_convective_source + use inifis_mod + use xios_output_mod, only: send_xios_field + implicit none + +! 0. DECLARATIONS: + + ! 0.1 INPUTS + + INTEGER, intent(in):: ngrid, nlayer + REAL, intent(in):: dtime ! Time step of the Physics + REAL, intent(in):: zmax_therm(ngrid) ! Altitude of max velocity thermals (m) + REAL, intent(in):: pp(ngrid, nlayer) ! Pressure at full levels + REAL, intent(in):: pt(ngrid, nlayer) ! Temperature at full levels + REAL, intent(in):: pu(ngrid, nlayer) ! Zonal wind at full levels + REAL, intent(in):: pv(ngrid, nlayer) ! Meridional wind at full levels + REAL, intent(in):: pdt(ngrid, nlayer) ! Tendency on temperature + REAL, intent(in):: pdu(ngrid, nlayer) ! Tendency on zonal wind + REAL, intent(in):: pdv(ngrid, nlayer) ! Tendency on meridional wind + + ! 0.2 OUTPUTS + + REAL, intent(out):: zustr(ngrid) ! Zonal surface stress + REAL, intent(out):: zvstr(ngrid) ! Meridional surface stress + REAL, intent(out):: d_t(ngrid, nlayer) ! Tendency on temperature + REAL, intent(out):: d_u(ngrid, nlayer) ! Tendency on zonal wind + REAL, intent(out):: d_v(ngrid, nlayer) ! Tendency on meridional wind + + ! 0.3 INTERNAL ARRAYS + + REAL :: tt(ngrid, nlayer) ! Temperature at full levels + REAL :: uu(ngrid, nlayer) ! Zonal wind at full levels + REAL :: vv(ngrid, nlayer) ! Meridional wind at full levels + REAL :: bvlow(ngrid) ! Qu est-ce ? + REAL :: dz + + INTEGER ii, jj, ll + + ! 0.3.0 Time scale of the like cycle of the waves parametrized + REAL, parameter:: deltat = 24. * 3600. + + ! 0.3.1 Gravity waves specifications + INTEGER, parameter:: nk = 2 ! number of horizontal wavenumbers + INTEGER, parameter:: np = 2 ! directions (eastward and westward) phase speed + INTEGER, parameter:: no = 2 ! absolute values of phase speed + INTEGER, parameter:: na = 5 ! Number of realizations to get the phase speed + INTEGER, parameter:: nw = nk * np *no ! Total numbers of gravity waves + INTEGER jk, jp, jo, jw + REAL, parameter:: kmax = 7.e-4 ! Max horizontal wavenumber + REAL, parameter:: kmin = 2.e-5 ! Min horizontal wavenumber + !REAL, parameter:: cmax = 30. ! Max horizontal absolute phase velocity + REAL, parameter:: cmax = 100. ! Test for Saturn: Max horizontal absolute phase velocity + REAL, parameter:: cmin = 1. ! Min horizontal absolute phase velocity + REAL cpha ! absolute phase velocity frequency + REAL zk(nw, ngrid) ! horizontal wavenumber amplitude + REAL zp(nw, ngrid) ! horizontal wavenumber angle + REAL zo(nw, ngrid) ! absolute frequency + + REAL intr_freq_m(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev below the full level (previous name: ZOM) + REAL intr_freq_p(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev above the full level (previous name: ZOP) + REAL wwm(nw, ngrid) ! Wave EP-fluxes at the 1/2 level below the full level + REAL wwp(nw, ngrid) ! Wave EP-fluxes at the 1/2 level above the full level + REAL u_epflux_p(nw, ngrid) ! Partial zonal flux (=for each wave) at the 1/2 level above the full level (previous name: RUWP) + REAL v_epflux_p(nw, ngrid) ! Partial meridional flux (=for each wave) at the 1/2 level above the full level (previous name: RVWP) + REAL u_epflux_tot(ngrid, nlayer + 1) ! Total zonal flux (=for all waves (nw)) at the 1/2 level above the full level (3D) (previous name: RUW) + REAL v_epflux_tot(ngrid, nlayer + 1) ! Total meridional flux (=for all waves (nw)) at the 1/2 level above the full level (3D) (previous name: RVW) + REAL epflux_0(nw, ngrid) ! Fluxes at launching level (previous name: RUW0) + INTEGER launch ! Launching altitude + REAL, parameter:: xlaunch = 0.4 ! Control the launching altitude + REAL, parameter:: zmaxth_top = 8000. ! Top of convective layer (approx.) + + REAL prec(ngrid) ! precipitations + REAL prmax ! Max value of precipitation + + ! 0.3.2 Parameters of waves dissipations + REAL, parameter:: sat = 1. ! saturation parameter + REAL, parameter:: rdiss = 1. ! coefficient of dissipation + REAL, parameter:: zoisec = 1.e-6 ! security for intrinsic freguency + + ! 0.3.3 Background flow at 1/2 levels and vertical coordinate + REAL H0bis(ngrid, nlayer) ! characteristic height of the atmosphere + REAL, save:: H0 ! characteristic height of the atmosphere + REAL, parameter:: pr = 250 ! Reference pressure [Pa] + REAL, parameter:: tr = 220. ! Reference temperature [K] + REAL zh(ngrid, nlayer + 1) ! Log-pressure altitude (constant H0) + REAL zhbis(ngrid, nlayer + 1) ! Log-pressure altitude (varying H) + REAL uh(ngrid, nlayer + 1) ! Zonal wind at 1/2 levels + REAL vh(ngrid, nlayer + 1) ! Meridional wind at 1/2 levels + REAL ph(ngrid, nlayer + 1) ! Pressure at 1/2 levels + REAL, parameter:: psec = 1.e-6 ! Security to avoid division by 0 pressure + REAL bv(ngrid, nlayer + 1) ! Brunt Vaisala freq. at 1/2 levels + REAL, parameter:: bvsec = 1.e-5 ! Security to avoid negative BV + REAL href(nlayer + 1) ! Reference altitude for launching altitude + + + + REAL ran_num_1, ran_num_2, ran_num_3 + + + LOGICAL, save :: firstcall = .true. + +!-------------------------------------------------------------------------------------------------- +! 1. INITIALISATION +!------------------ + IF (firstcall) THEN + write(*,*) "nonoro_gwd_ran: FLott non-oro GW scheme is active!" + ! Characteristic vertical scale height + H0 = r * tr / g + ! Control + if (deltat .LT. dtime) THEN + call abort_physic("nonoro_gwd_ran","gwd random: deltat lower than dtime!",1) + endif + if (nlayer .LT. nw) THEN + call abort_physic("nonoro_gwd_ran","gwd random: nlayer lower than nw!",1) + endif + firstcall = .false. + ENDIF + gwd_convective_source = .false. + + ! Compute subroutine's current values of temperature and winds + tt(:,:) = pt(:,:) + dtime * pdt(:,:) + uu(:,:) = pu(:,:) + dtime * pdu(:,:) + vv(:,:) = pv(:,:) + dtime * pdv(:,:) +! print*,'epflux_max just after firstcall:', epflux_max + + +! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS +!---------------------------------------------------- + ! Pressure and Inv of pressure, temperature at 1/2 level + DO ll = 2, nlayer + ph(:, ll) = exp((log(pp(:, ll)) + log(pp(:, ll - 1))) / 2.) + end DO + + ph(:, nlayer + 1) = 0. + ph(:, 1) = 2. * pp(:, 1) - ph(:, 2) + + ! Launching altitude for reproductible case + DO ll = 2, nlayer + href(ll) = exp((log(presnivs(ll)) + log(presnivs(ll - 1))) / 2.) + end DO + href(nlayer + 1) = 0. + href(1) = 2. * presnivs(1) - href(2) + + launch = 0. + DO ll =1, nlayer + IF (href (ll) / href(1) > xlaunch) launch = ll + end DO + + ! Log-pressure vertical coordinate + DO ll = 1, nlayer + 1 + zh(:,ll) = H0 * log(pr / (ph(:,ll) + psec)) + end DO + + ! Winds and Brunt Vaisala frequency + DO ll = 2, nlayer + uh(:, ll) = 0.5 * (uu(:, ll) + uu(:, ll - 1)) ! Zonal wind + vh(:, ll) = 0.5 * (vv(:, ll) + vv(:, ll - 1)) ! Meridional wind + + bv(:, ll) = 0.5 * (tt(:, ll) + tt(:, ll - 1)) & + * r**2 / cpp / H0**2 & + + (tt(:, ll) - tt(:, ll - 1)) & + / (zh(:, ll) - zh(:, ll - 1)) & + * r / H0 + bv(:, ll) = sqrt(max(bvsec, bv(:,ll))) + end DO + + bv(:, 1) = bv(:, 2) + uh(:, 1) = 0. + vh(:, 1) = 0. + bv(:, nlayer + 1) = bv(:, nlayer) + uh(:, nlayer + 1) = uu(:, nlayer) + vh(:, nlayer + 1) = vv(:, nlayer) + +! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY +!----------------------------------------- +! The mod function of here a weird arguments +! are used to produce the waves characteristics +! in a stochastic way + DO jw = 1, nw + DO ii = 1, ngrid + + ran_num_1 = mod(tt(ii, jw) * 10., 1.) + ran_num_2 = mod(tt(ii, jw) * 100., 1.) + + ! angle (0 or pi so far) + zp(jw, ii) = (sign(1., 0.5 - ran_num_1) & + + 1.) * pi / 2. + ! horizontal wavenumber amplitude + zk(jw, ii) = kmin + (kmax - kmin) * ran_num_2 + ! horizontal phase speed + cpha = 0. + DO jj = 1, na + ran_num_3 = mod(tt(ii, jw + 3 * jj)**2, 1.) + cpha = cpha + 2. * cmax * & + (ran_num_3 - 0.5) * & + sqrt(3.) / sqrt(na * 1.) + end DO + IF (cpha < 0.) THEN + cpha = - 1. * cpha + zp (jw, ii) = zp(jw, ii) + pi + ENDIF + ! Intrinsic frequency + zo(jw, ii) = cpha * zk(jw, ii) + ! Intrinsic frequency is imposed + 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 level + epflux_0(jw, ii) = epflux_max & + * mod(100. * (uu(ii, jw)**2 + vv(ii, jw)**2), 1.) + + end DO + end DO + +! print*,'epflux_max just after waves charac. ramdon:', epflux_max +! print*,'epflux_0 just after waves charac. ramdon:', epflux_0 +! 4. COMPUTE THE FLUXES +!---------------------- + ! 4.1 Vertical velocity at launching altitude to ensure + ! the correct value to the imposed fluxes. + !------------------------------------------------------ + DO jw = 1, nw + ! Evaluate intrinsic frequency at launching altutide: + intr_freq_p(jw, :) = zo(jw, :) & + - zk(jw, :) * cos(zp(jw, :)) * uh(:, launch) & + - zk(jw, :) * sin(zp(jw, :)) * vh(:, launch) + end DO + + IF (gwd_convective_source) THEN + DO jw = 1, nw + ! VERSION WITH CONVECTIVE SOURCE ! designed for Earth + + ! Vertical velocity at launch level, value to ensure the + ! imposed mmt flux factor related to the convective forcing: + ! precipitations. + + ! tanh limitation to values above prmax: +! WWP(JW, :) = epflux_0(JW, :) & +! * (r / cpp / H0 * RLVTT * PRMAX * TANH(PREC(:) / PRMAX))**2 +! Here, we neglected the kinetic energy providing of the thermodynamic +! phase change + +! + + ! Factor related to the characteristics of the waves: + wwp(jw, :) = wwp(jw, :) * zk(jw, :)**3 / kmin / bvlow(:) & + / MAX(ABS(intr_freq_p(jw, :)), zoisec)**3 + + ! Moderation by the depth of the source (dz here): + wwp(jw, :) = wwp(jw, :) & + * exp(- bvlow(:)**2 / max(abs(intr_freq_p(jw, :)), zoisec)**2 & + * zk(jw, :)**2 * dz**2) + + ! Put the stress in the right direction: + u_epflux_p(jw, :) = intr_freq_p(jw, :) / max(abs(intr_freq_p(jw, :)), zoisec)**2 & + * bv(:, launch) * cos(zp(jw, :)) * wwp(jw, :)**2 + v_epflux_p(JW, :) = intr_freq_p(jw, :) / max(abs(intr_freq_p(jw, :)), zoisec)**2 & + * bv(:, launch) * sin(zp(jw, :)) * wwp(jw, :)**2 + end DO + ELSE ! VERSION WITHOUT CONVECTIVE SOURCE + ! Vertical velocity at launch level, value to ensure the imposed + ! mom flux: + DO jw = 1, nw + ! WW is directly a flux, here, not vertical velocity anymore + wwp(jw, :) = epflux_0(JW,:) + u_epflux_p(jw, :) = cos(zp(jw, :)) * sign(1., intr_freq_p(jw, :)) * epflux_0(jw, :) + v_epflux_p(jw, :) = sin(zp(jw, :)) * sign(1., intr_freq_p(jw, :)) * epflux_0(jw, :) + + end DO + ENDIF + + ! 4.2 Initial flux at launching altitude + !--------------------------------------- + u_epflux_tot(:, launch) = 0. + v_epflux_tot(:, launch) = 0. + DO jw = 1, nw + u_epflux_tot(:, launch) = u_epflux_tot(:, launch) + u_epflux_p(jw, :) + v_epflux_tot(:, launch) = v_epflux_tot(:, launch) + v_epflux_p(jw, :) + end DO + + ! 4.3 Loop over altitudes, with passage from one level to the next done by: + !-------------------------------------------------------------------------- + ! i) conserving the EP flux, + ! ii) dissipating a little, + ! iii) testing critical levels, + ! iv) testing the breaking. + !-------------------------------------------------------------------------- + DO ll = launch, nlayer - 1 + ! Warning! all the physics is here (passage from one level to the next + DO jw = 1, nw + intr_freq_m(jw, :) = intr_freq_p(jw, :) + wwm(jw, :) = wwp(jw, :) + ! Intrinsic frequency + intr_freq_p(jw, :) = zo(jw, :) - zk(jw, :) * cos(zp(jw, :)) * uh(:, ll + 1) & + - zk(jw, :) * sin(zp(jw, :)) * vh(:, ll + 1) + ! Vertical velocity in flux formulation + wwp(jw, :) = min( & + ! No breaking (Eq. 6): + wwm(jw, :) & + ! Dissipation (Eq. 8): + * exp(-rdiss * pr / (ph(:, ll + 1) + ph (:, ll)) & + * ((bv(:, ll + 1) + bv (:, ll)) / 2.)**3 & + / max(abs(intr_freq_p(jw, :) + intr_freq_m(jw, :)) / 2., zoisec)**4 & + * zk(jw, :)**3 * (zh(:, ll + 1) - zh(:, ll))) , & + ! Critical levels (forced to zero if intrinsic frequency + ! changes sign) + max(0., sign(1., intr_freq_p(jw, :) * intr_freq_m(jw, :))) & + ! Saturation (Eq. 12) + * abs(intr_freq_p(jw, :))**3 / bv(:, ll + 1) & + * exp(-zh(:, ll + 1) / H0) & + * sat**2 * kmin**2 / zk(jw, :)**4) + end DO + + ! Evaluate EP-flux from Eq. 7 and give the right orientation to the stress + DO jw = 1, nw + u_epflux_p(jw, :) = sign(1., intr_freq_p(jw, :)) * cos(zp(jw, :)) * wwp(jw, :) + v_epflux_p(jw, :) = sign(1., intr_freq_p(jw, :)) * sin(zp(jw, :)) * wwp(jw, :) + end DO + + u_epflux_tot(:, ll + 1) = 0. + v_epflux_tot(:, ll + 1) = 0. + DO jw = 1, nw + u_epflux_tot(:, ll + 1) = u_epflux_tot(:, ll + 1) + u_epflux_p(jw, :) + v_epflux_tot(:, ll + 1) = v_epflux_tot(:, ll + 1) + v_epflux_p(jw, :) + east_gwstress(:, ll) = east_gwstress(:, ll) & + + max(0., u_epflux_p(jw, :)) / float(nw) + west_gwstress(:, ll) = west_gwstress(:, ll) & + + min(0., u_epflux_p(jw, ll)) / float(nw) + end DO + end DO ! DO LL = LAUNCH, nlayer - 1 + +! print*,'u_epflux_tot just after launching:', u_epflux_tot +! print*,'v_epflux_tot just after launching:', v_epflux_tot +! print*,'u_epflux_p just after launching:', maxval(u_epflux_p), minval(u_epflux_p) +! print*,'v_epflux_p just after launching:', maxval(v_epflux_p), minval(v_epflux_p) + +! 5. TENDENCY CALCULATIONS +!------------------------- + ! 5.1 Flow rectification at the top and in the low layers + ! -------------------------------------------------------- + ! Warning, this is the total on all GW... + + u_epflux_tot(:, nlayer + 1) = 0. + v_epflux_tot(:, nlayer + 1) = 0. + + ! Here, big change compared to FLott version: + ! We compensate (u_epflux_(:, LAUNCH), ie total emitted upward flux + ! over the layers max(1,LAUNCH-3) to LAUNCH-1 + DO LL = 1, max(1,LAUNCH-3) + u_epflux_tot(:, LL) = 0. + v_epflux_tot(:, LL) = 0. + end DO + DO LL = max(2,LAUNCH-2), LAUNCH-1 + u_epflux_tot(:, LL) = u_epflux_tot(:, LL - 1) & + + u_epflux_tot(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) & + / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) + v_epflux_tot(:, LL) = v_epflux_tot(:, LL - 1) & + + v_epflux_tot(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) & + / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) + east_gwstress(:,LL) = east_gwstress(:, LL - 1) & + + east_gwstress(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) & + / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) + west_gwstress(:,LL) = west_gwstress(:, LL - 1) & + + west_gwstress(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) & + / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) + end DO + ! This way, the total flux from GW is zero, but there is a net transport + ! (upward) that should be compensated by circulation + ! and induce additional friction at the surface + + ! 5.2 AR-1 RECURSIVE FORMULA (13) IN VERSION 4 + !--------------------------------------------- + DO LL = 1, nlayer + d_u(:, LL) = (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL))! & + !d_u(:, LL) = G * (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL))! & + ! / (PH(:, LL + 1) - PH(:, LL)) * DTIME + d_v(:, LL) = (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL))! & + !d_v(:, LL) = G * (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL))! & + ! / (PH(:, LL + 1) - PH(:, LL)) * DTIME + ENDDO + d_t(:,:) = 0. + + ! 5.3 Update tendency of wind with the previous (and saved) values + !----------------------------------------------------------------- + d_u(:,:) = DTIME/DELTAT/REAL(NW) * d_u(:,:) & + + (1.-DTIME/DELTAT) * du_nonoro_gwd(:,:) + d_v(:,:) = DTIME/DELTAT/REAL(NW) * d_v(:,:) & + + (1.-DTIME/DELTAT) * dv_nonoro_gwd(:,:) + du_nonoro_gwd(:,:) = d_u(:,:) + dv_nonoro_gwd(:,:) = d_v(:,:) + +! print*,'u_epflux_tot just after tendency:', u_epflux_tot +! print*,'v_epflux_tot just after tendency:', v_epflux_tot +! print*,'d_u just after tendency:', maxval(d_u(:,:)), minval(d_u) +! print*,'d_v just after tendency:', maxval(d_v(:,:)), minval(d_v) + ! Cosmetic: evaluation of the cumulated stress + + ZUSTR(:) = 0. + ZVSTR(:) = 0. + DO LL = 1, nlayer + ZUSTR(:) = ZUSTR(:) + D_U(:, LL) !/ g * (PH(:, LL + 1) - PH(:, LL)) + ZVSTR(:) = ZVSTR(:) + D_V(:, LL) !/ g * (PH(:, LL + 1) - PH(:, LL)) + ENDDO + + call send_xios_field("du_nonoro", d_u) + call send_xios_field("dv_nonoro", d_v) + call send_xios_field("east_gwstress", east_gwstress) + call send_xios_field("west_gwstress", west_gwstress) + + + END SUBROUTINE nonoro_gwd_ran + +! =================================================================== +! Subroutines used to write variables of memory in start files +! =================================================================== + + SUBROUTINE ini_nonoro_gwd_ran(ngrid,nlayer) + + IMPLICIT NONE + + INTEGER, INTENT (in) :: ngrid ! number of atmospheric columns + INTEGER, INTENT (in) :: nlayer ! number of atmospheric layers + + allocate(du_nonoro_gwd(ngrid, nlayer)) + allocate(dv_nonoro_gwd(ngrid, nlayer)) + allocate(east_gwstress(ngrid, nlayer)) + allocate(west_gwstress(ngrid, nlayer)) + + END SUBROUTINE ini_nonoro_gwd_ran +! ---------------------------------- + SUBROUTINE end_nonoro_gwd_ran + + IMPLICIT NONE + + if (allocated(du_nonoro_gwd)) deallocate(du_nonoro_gwd) + if (allocated(dv_nonoro_gwd)) deallocate(dv_nonoro_gwd) + if (allocated(east_gwstress)) deallocate(east_gwstress) + if (allocated(west_gwstress)) deallocate(west_gwstress) + + END SUBROUTINE end_nonoro_gwd_ran + +END MODULE nonoro_gwd_ran_mod + Index: /trunk/LMDZ.GENERIC/libf/phystd/phyetat0_mod.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/phyetat0_mod.F90 (revision 2298) +++ /trunk/LMDZ.GENERIC/libf/phystd/phyetat0_mod.F90 (revision 2299) @@ -22,5 +22,6 @@ use slab_ice_h, only: noceanmx use callkeys_mod, only: CLFvarying,surfalbedo,surfemis - + use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd, & + east_gwstress, west_gwstress implicit none @@ -389,4 +390,50 @@ call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime) endif ! of if (startphy_file) + +! Non-orographic gravity waves +if (startphy_file) then + call get_field("du_nonoro_gwd",du_nonoro_gwd,found,indextime) + if (.not.found) then + write(*,*) "phyetat0: not in file" + du_nonoro_gwd(:,:)=0. + else + write(*,*) "phyetat0: Memory of zonal wind tendency due to non-orographic GW" + write(*,*) " range:", & + minval(du_nonoro_gwd), maxval(du_nonoro_gwd) + endif +endif ! of if (startphy_file) +if (startphy_file) then + call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found,indextime) + if (.not.found) then + write(*,*) "phyetat0: not in file" + dv_nonoro_gwd(:,:)=0. + else + write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW" + write(*,*) " range:", & + minval(dv_nonoro_gwd), maxval(dv_nonoro_gwd) + endif +endif ! of if (startphy_file) +if (startphy_file) then + call get_field("east_gwstress",east_gwstress,found,indextime) + if (.not.found) then + write(*,*) "phyetat0: not in file" + east_gwstress(:,:)=0. + else + write(*,*) "phyetat0: Memory of eastward stress profile due to non-orographic GW" + write(*,*) " range:", & + minval(east_gwstress), maxval(east_gwstress) + endif +endif ! of if (startphy_file) +if (startphy_file) then + call get_field("west_gwstress",west_gwstress,found,indextime) + if (.not.found) then + write(*,*) "phyetat0: not in file" + west_gwstress(:,:)=0. + else + write(*,*) "phyetat0: Memory of westward stress profile due to non-orographic GW" + write(*,*) " range:", & + minval(west_gwstress), maxval(west_gwstress) + endif +endif ! of if (startphy_file) ! ! close file: Index: /trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90 (revision 2298) +++ /trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90 (revision 2299) @@ -142,5 +142,7 @@ use tracer_h, only: noms use slab_ice_h, only: noceanmx - use callkeys_mod, only: ok_slab_ocean + use callkeys_mod, only: ok_slab_ocean, calllott_nonoro + use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd, & + east_gwstress, west_gwstress implicit none @@ -208,4 +210,13 @@ enddo endif ! of if (nq>0) + +! Non-orographic gavity waves +if (calllott_nonoro) then + call put_field("du_nonoro_gwd","Zonal wind tendency due to GW",du_nonoro_gwd) + call put_field("dv_nonoro_gwd","Meridional wind tendency due to GW",dv_nonoro_gwd) + call put_field("east_gwstress","Eastward stress profile due to GW",east_gwstress) + call put_field("west_gwstress","Westward stress profile due to GW",west_gwstress) +endif + ! close file CALL close_restartphy Index: /trunk/LMDZ.GENERIC/libf/phystd/phys_state_var_mod.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/phys_state_var_mod.F90 (revision 2298) +++ /trunk/LMDZ.GENERIC/libf/phystd/phys_state_var_mod.F90 (revision 2299) @@ -18,4 +18,5 @@ zmea, zstd, zsig, zgam, zthe use turb_mod, only: q2,sensibFlux,wstar,ustar,tstar,hfmax_th,zmax_th + use nonoro_gwd_ran_mod, only: ini_nonoro_gwd_ran, end_nonoro_gwd_ran real,allocatable,dimension(:,:),save :: ztprevious ! Previous loop Atmospheric Temperature (K) @@ -179,5 +180,6 @@ allocate(hfmax_th(klon)) allocate(zmax_th(klon)) - + ! allocate arrays in "nonoro_gwd_ran_mod" + call ini_nonoro_gwd_ran(klon,klev) END SUBROUTINE phys_state_var_init @@ -251,5 +253,6 @@ deallocate(hfmax_th) deallocate(zmax_th) - + ! deallocate arrays in "nonoro_gwd_ran_mod" + call end_nonoro_gwd_ran END SUBROUTINE phys_state_var_end Index: /trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90 (revision 2298) +++ /trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90 (revision 2299) @@ -48,4 +48,5 @@ use time_phylmdz_mod, only: daysec use callkeys_mod + use nonoro_gwd_ran_mod, only: nonoro_gwd_ran use conc_mod use phys_state_var_mod @@ -410,4 +411,9 @@ integer num_run logical,save :: ice_update +! Non-oro GW tendencies + REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer) + REAL d_t_hin(ngrid,nlayer) +! Diagnostics 2D of gw_nonoro + REAL zustrhi(ngrid), zvstrhi(ngrid) @@ -497,11 +503,10 @@ endif #endif - if (pday.ne.day_ini) then - write(*,*) "ERROR in physiq.F90:" - write(*,*) "bad synchronization between physics and dynamics" - write(*,*) "dynamics day: ",pday - write(*,*) "physics day: ",day_ini - stop + write(*,*) "ERROR in physiq.F90:" + write(*,*) "bad synchronization between physics and dynamics" + write(*,*) "dynamics day: ",pday + write(*,*) "physics day: ",day_ini + stop endif @@ -680,6 +685,6 @@ if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d. call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, & - ptimestep,pday+nday,time_phys,cell_area, & - albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe) + ptimestep,pday+nday,time_phys,cell_area, & + albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe) endif #endif @@ -1239,4 +1244,32 @@ endif ! end of 'calladj' +!---------------------------------------------- +! Non orographic Gravity Waves: +!--------------------------------------------- + IF (calllott_nonoro) THEN + + CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep,pplay, & + zmax_th, &! max altitude reached by thermals (m) + pt, pu, pv, & + pdt, pdu, pdv, & + zustrhi,zvstrhi, & + d_t_hin, d_u_hin, d_v_hin) + +! Update tendencies + pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) & + + d_t_hin(1:ngrid,1:nlayer) +! d_t_hin(:,:)= d_t_hin(:,:)/ptimestep ! K/s (for outputs?) + pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) & + + d_u_hin(1:ngrid,1:nlayer) +! d_u_hin(:,:)= d_u_hin(:,:)/ptimestep ! (m/s)/s (for outputs?) + pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) & + + d_v_hin(1:ngrid,1:nlayer) +! d_v_hin(:,:)= d_v_hin(:,:)/ptimestep ! (m/s)/s (for outputs?) + print*,'du_nonoro: max ', maxval(d_u_hin), 'min ', minval(d_u_hin) + print*,'dv_nonoro: max ', maxval(d_v_hin), 'min ', minval(d_v_hin) + + ENDIF ! of IF (calllott_nonoro) + + !----------------------------------------------- @@ -2112,4 +2145,10 @@ call writediagfi(ngrid,'fraca','Updraft fraction','', 3, fraca) endif + +! GW non-oro outputs + if (calllott_nonoro) then + call WRITEDIAGFI(ngrid,"dugwno","GW non-oro dU","m/s2", 3,d_u_hin/ptimestep) + call WRITEDIAGFI(ngrid,"dvgwno","GW non-oro dV","m/s2", 3,d_v_hin/ptimestep) + endif ! Total energy balance diagnostics @@ -2189,5 +2228,5 @@ !call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw) !call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw) - !call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad) + call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad) ! call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)