Ignore:
Timestamp:
Dec 16, 2021, 3:09:04 PM (3 years ago)
Author:
emillour
Message:

Generic GCM:
Fixes and improvements in the Non-orographic GW scheme, namely:

  • increments are not tendencies
  • missing rho factor in EP flux computation
  • missing rho at launch altitude
  • changed inputs, because R and Cp are needed to compute rho and BV

JL

Location:
trunk/LMDZ.GENERIC/libf/phystd
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90

    r2542 r2595  
    6464!$OMP THREADPRIVATE(photochem,photoheat,jonline,depos)
    6565      logical,save :: calllott_nonoro
    66       logical,save :: gwd_convective_source
    67 !$OMP THREADPRIVATE(calllott_nonoro,gwd_convective_source)
     66!$OMP THREADPRIVATE(calllott_nonoro)
    6867      logical,save :: global1d
    6968      real,save    :: szangle
    … …  
    145144      real,save :: surfemis
    146145!$OMP THREADPRIVATE(surfalbedo,surfemis)
    147       real,save :: epflux_max
    148 !$OMP THREADPRIVATE(epflux_max)
    149146      real,save :: noseason_day
    150147!$OMP THREADPRIVATE(noseason_day)

  • trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90

    r2583 r2595  
    329329       ": calllott_nonoro = ",calllott_nonoro
    330330     
    331      if (is_master) write(*,*)trim(rname)//&
    332        ": Max Value of EP flux at launching altitude"
    333      epflux_max=0.0                ! default value
    334      call getin_p("epflux_max",epflux_max)
    335      if (is_master) write(*,*)trim(rname)//": epflux_max = ",epflux_max
    336 
    337331     if (is_master) write(*,*) trim(rname)//": call thermal plume model ?"
    338332     calltherm=.false. ! default value

  • trunk/LMDZ.GENERIC/libf/phystd/nonoro_gwd_ran_mod.F90

    r2550 r2595  
    1010CONTAINS
    1111
    12       SUBROUTINE nonoro_gwd_ran(ngrid, nlayer, dtime, pp, &
     12      SUBROUTINE nonoro_gwd_ran(ngrid, nlayer, dtime, cpnew, rnew, pp, &
    1313                   zmax_therm, pt, pu, pv, pdt, pdu, pdv, &
    1414                   zustr, zvstr, d_t, d_u, d_v)
    … …  
    3838!
    3939! ADAPTED FOR GENERIC  D.BARDET    01/2020
     40! UPDATED              J.LIU       12/2021  The rho (density) at the specific
     41!                                           locations is introduced. The equation
     42!                                           of EP-flux is corrected by adding the
     43!                                           term of density at launch (source)
     44!                                           altitude(level). Bugs in BV,d_u,d_v fixed.
    4045!---------------------------------------------------------------------------------
    4146
    4247
    4348
    44      use comcstfi_mod, only: g, pi, cpp, r
     49     use comcstfi_mod, only: g, pi, r
    4550     USE ioipsl_getin_p_mod, ONLY : getin_p
    46      use assert_m, only : assert
    4751     use vertical_layers_mod, only : presnivs
    48      use callkeys_mod, only : epflux_max, & ! Max EP flux value at launching altitude (previous name: RUWMAX)
    49                               gwd_convective_source
    50      use inifis_mod
    5152     use geometry_mod, only: cell_area
    5253#ifdef CPP_XIOS
    5354     use xios_output_mod, only: send_xios_field
    5455#endif
    55  implicit none
    56 
    57 ! 0. DECLARATIONS:
     56     implicit none
     57
     58     ! 0. DECLARATIONS:
    5859
    5960     ! 0.1 INPUTS
    6061
    61      INTEGER, intent(in):: ngrid, nlayer
     62     INTEGER, intent(in):: ngrid           ! number of atmospheric columns
     63     INTEGER, intent(in):: nlayer          ! number of atmospheric columns
    6264     REAL, intent(in):: dtime              ! Time step of the Physics
    6365     REAL, intent(in):: zmax_therm(ngrid)  ! Altitude of max velocity thermals (m)
    64      REAL, intent(in):: pp(ngrid, nlayer)  ! Pressure at full levels
    65      REAL, intent(in):: pt(ngrid, nlayer)  ! Temperature at full levels
    66      REAL, intent(in):: pu(ngrid, nlayer)  ! Zonal wind at full levels
    67      REAL, intent(in):: pv(ngrid, nlayer)  ! Meridional wind at full levels
    68      REAL, intent(in):: pdt(ngrid, nlayer) ! Tendency on temperature
    69      REAL, intent(in):: pdu(ngrid, nlayer) ! Tendency on zonal wind
    70      REAL, intent(in):: pdv(ngrid, nlayer) ! Tendency on meridional wind
     66     REAL,INTENT(IN) :: cpnew(ngrid,nlayer)! Cp of the atmosphere
     67     REAL,INTENT(IN) :: rnew(ngrid,nlayer) ! R of the atmosphere
     68     REAL, intent(in):: pp(ngrid, nlayer)  ! Pressure at full levels(Pa)
     69     REAL, intent(in):: pt(ngrid, nlayer)  ! Temperature at full levels(K)
     70     REAL, intent(in):: pu(ngrid, nlayer)  ! Zonal wind at full levels(m/s)
     71     REAL, intent(in):: pv(ngrid, nlayer)  ! Meridional wind at full levels(m/s)
     72     REAL, intent(in):: pdt(ngrid, nlayer) ! Tendency on temperature(K/s)
     73     REAL, intent(in):: pdu(ngrid, nlayer) ! Tendency on zonal wind(m/s/s)
     74     REAL, intent(in):: pdv(ngrid, nlayer) ! Tendency on meridional wind(m/s/s)
    7175     
    7276     ! 0.2 OUTPUTS
    … …  
    7478     REAL, intent(out):: zustr(ngrid)         ! Zonal surface stress
    7579     REAL, intent(out):: zvstr(ngrid)         ! Meridional surface stress
    76      REAL, intent(out):: d_t(ngrid, nlayer)   ! Tendency on temperature
    77      REAL, intent(out):: d_u(ngrid, nlayer)   ! Tendency on zonal wind
    78      REAL, intent(out):: d_v(ngrid, nlayer)   ! Tendency on meridional wind
     80     REAL, intent(out):: d_t(ngrid, nlayer)   ! Tendency on temperature (K/s) due to gravity waves (not used set to zero)
     81     REAL, intent(out):: d_u(ngrid, nlayer)   ! Tendency on zonal wind (m/s/s) due to gravity waves
     82     REAL, intent(out):: d_v(ngrid, nlayer)   ! Tendency on meridional wind (m/s/s) due to gravity waves
    7983
    8084     ! 0.3 INTERNAL ARRAYS
    8185
    8286     REAL :: tt(ngrid, nlayer) ! Temperature at full levels
     87     REAL :: rho(ngrid, nlayer)  ! Mass density at full levels
    8388     REAL :: uu(ngrid, nlayer) ! Zonal wind at full levels
    8489     REAL :: vv(ngrid, nlayer) ! Meridional wind at full levels
    85      REAL :: bvlow(ngrid)      ! Qu est-ce ?
    86      REAL :: dz
     90     REAL :: bvlow(ngrid)      ! initial N at each grid (not used)
     91     REAL :: dz                ! depth of the GW source if gwd_convective_source=.true.
    8792
    8893     INTEGER ii, jj, ll
    … …  
    98103     INTEGER, parameter:: nw = nk * np *no ! Total numbers of gravity waves
    99104     INTEGER jk, jp, jo, jw
    100      REAL kstar                  ! Control value to constrain the min horizontal wavenumber by the horizontal resolution
    101      REAL, parameter:: kmax = 1.e-4 ! Max horizontal wavenumber=N/u, u(=30) zonal wind velocity at launch
    102      !REAL, parameter:: kmax = 4.e-5 ! Max horizontal wavenumber=N/u, u(=70) zonal wind velocity at launch
    103      REAL, parameter:: kmin = 2.e-6 ! Min horizontal wavenumber=1/sqrt(DxDy) Dx and Dy horizontal grid
    104      REAL, parameter:: cmax = 30.   ! Max horizontal absolute phase velocity=zonal wind at launch
     105
     106     REAL kstar                      ! Control value to constrain the min horizontal wavenumber by the horizontal resolution
     107     REAL, save :: kmax ! Max horizontal wavenumber=N/u, u(=30) zonal wind velocity at launch,lambda min,lambda=2*pi/kmax=62.8 km
     108!$OMP THREADPRIVATE(kmax)
     109     !REAL, parameter:: kmax = 4.e-5 ! Max horizontal wavenumber=N/u, u(=70) zonal wind velocity at launch,lambda min,lambda=2*pi/kmax
     110     REAL, parameter:: kmin = 2.e-6  ! Min horizontal wavenumber=1/sqrt(DxDy) Dx and Dy horizontal grid,lambda max,lambda=2*pi/kmax=314 km
     111     REAL, save :: cmax     ! Max horizontal absolute phase velocity=zonal wind at launch
     112!$OMP THREADPRIVATE(cmax)
    105113     !REAL, parameter:: cmax = 100.  ! Test for Saturn: Max horizontal absolute phase velocity
    106114     !REAL, parameter:: cmax = 70.   ! Test for Saturn: Max horizontal absolute phase velocity
    107      REAL, parameter:: cmin = 1.    ! Min horizontal absolute phase velocity
    108      REAL cpha                      ! absolute phase velocity frequency
    109      REAL zk(nw, ngrid)             ! horizontal wavenumber amplitude
    110      REAL zp(nw, ngrid)             ! horizontal wavenumber angle
    111      REAL zo(nw, ngrid)             ! absolute frequency
     115     REAL, parameter:: cmin = 1.     ! Min horizontal absolute phase velocity(not used)
     116     REAL cpha                       ! absolute phase velocity frequency
     117     REAL zk(nw, ngrid)              ! horizontal wavenumber amplitude
     118     REAL zp(nw, ngrid)              ! horizontal wavenumber angle
     119     REAL zo(nw, ngrid)              ! absolute frequency
    112120
    113121     REAL intr_freq_m(nw, ngrid)          ! Waves Intr. freq. at the 1/2 lev below the full level (previous name: ZOM)
    … …  
    120128     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)
    121129     REAL epflux_0(nw, ngrid)             ! Fluxes at launching level (previous name: RUW0)
     130     REAL, save :: epflux_max             ! Max EP flux value at launching altitude (previous name: RUWMAX, tunable)
     131!$OMP THREADPRIVATE(epflux_max)
    122132     INTEGER launch                       ! Launching altitude
    123      REAL, parameter:: xlaunch = 0.2      ! Control the launching altitude
     133     REAL, parameter:: xlaunch = 0.2      ! Control the launching altitude by pressure
    124134     REAL, parameter:: zmaxth_top = 8000. ! Top of convective layer in m (approx.)
    125135
    … …  
    128138
    129139     ! 0.3.2 Parameters of waves dissipations
    130      REAL, parameter:: sat   = 1.     ! saturation parameter
    131      REAL, parameter:: rdiss = 1.     ! coefficient of dissipation
    132      REAL, parameter:: zoisec = 1.e-6 ! security for intrinsic freguency
     140     REAL, save :: sat        ! saturation parameter(tunable)
     141!$OMP THREADPRIVATE(sat)
     142     REAL, save :: rdiss      ! coefficient of dissipation(tunable)
     143!$OMP THREADPRIVATE(rdiss)
     144     REAL, parameter:: zoisec = 1.e-10 ! security for intrinsic freguency
    133145
    134146     ! 0.3.3 Background flow at 1/2 levels and vertical coordinate
    135147     REAL H0bis(ngrid, nlayer)       ! characteristic height of the atmosphere
    136148     REAL, save::  H0                ! characteristic height of the atmosphere
     149!$OMP THREADPRIVATE(H0)
    137150     !REAL, parameter:: pr = 250      ! Reference pressure [Pa]
    138151     !REAL, parameter:: tr = 220.     ! Reference temperature [K]
    … …  
    145158     REAL vh(ngrid, nlayer + 1)      ! Meridional wind at 1/2 levels
    146159     REAL ph(ngrid, nlayer + 1)      ! Pressure at 1/2 levels
    147      REAL, parameter:: psec = 1.e-6  ! Security to avoid division by 0 pressure
     160     REAL, parameter:: psec = 1.e-12 ! Security to avoid division by 0 pressure(!!IMPORTANT: should be lower than the topmost layer's pressure)
    148161     REAL bv(ngrid, nlayer + 1)      ! Brunt Vaisala freq. at 1/2 levels
    149      REAL, parameter:: bvsec = 1.e-5 ! Security to avoid negative BV
     162     REAL, parameter:: bvsec = 1.e-6 ! Security to avoid negative BV (!!IMPORTANT: tunable, depends on which Planet)
    150163     REAL href(nlayer + 1)           ! Reference altitude for launching altitude
    151164
    … …  
    156169
    157170     LOGICAL, save :: firstcall = .true.
     171!$OMP THREADPRIVATE(firstcall)
     172     LOGICAL, save :: gwd_convective_source = .false.
     173!$OMP THREADPRIVATE(gwd_convective_source)
    158174
    159175!-------------------------------------------------------------------------------------------------- 
    160176! 1. INITIALISATION
    161 !------------------
     177!--------------------------------------------------------------------------------------------------
    162178     IF (firstcall) THEN
    163179        write(*,*) "nonoro_gwd_ran: FLott non-oro GW scheme is active!"
     180        epflux_max = 0.0 ! Default value !!
     181        call getin_p("nonoro_gwd_epflux_max", epflux_max)
     182        write(*,*) "nonoro_gwd_ran: epflux_max=", epflux_max
     183        sat = 1. ! default gravity waves saturation value !!
     184        call getin_p("nonoro_gwd_sat", sat)
     185        write(*,*) "nonoro_gwd_ran: sat=", sat     
     186        cmax = 30. ! default gravity waves phase velocity value !!
     187        call getin_p("nonoro_gwd_cmax", cmax)
     188        write(*,*) "nonoro_gwd_ran: cmax=", cmax
     189        rdiss=1 ! default coefficient of dissipation !!
     190        call getin_p("nonoro_gwd_rdiss", rdiss)
     191        write(*,*) "nonoro_gwd_ran: rdiss=", rdiss
     192        kmax = 1.e-4  ! default Max horizontal wavenumber !!
     193        call getin_p("nonoro_gwd_kmax", kmax)
     194        write(*,*) "nonoro_gwd_ran: kmax=", kmax
     195
    164196        ! Characteristic vertical scale height
    165197        H0 = r * tr / g
    … …  
    175207     gwd_convective_source = .false.
    176208
    177      ! Compute subroutine's current values of temperature and winds
     209! Compute subroutine's current values of temperature and winds
    178210     tt(:,:) = pt(:,:) + dtime * pdt(:,:)
    179211     uu(:,:) = pu(:,:) + dtime * pdu(:,:)
    180212     vv(:,:) = pv(:,:) + dtime * pdv(:,:)
     213! Compute the real mass density by rho=p/R(T)T
     214     DO ll=1,nlayer
     215       DO ii=1,ngrid
     216            rho(ii,ll) = pp(ii,ll)/(rnew(ii,ll)*tt(ii,ll))
     217       ENDDO
     218     ENDDO
    181219!    print*,'epflux_max just after firstcall:', epflux_max
    182220
    183 
     221!--------------------------------------------------------------------------------------------------
    184222! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
    185 !----------------------------------------------------
     223!--------------------------------------------------------------------------------------------------
    186224     ! Pressure and Inv of pressure, temperature at 1/2 level
    187225     DO ll = 2, nlayer
    … …  
    205243   
    206244     ! Log-pressure vertical coordinate
    207      DO ll = 1, nlayer + 1
    208         zh(:,ll) = H0 * log(pr / (ph(:,ll) + psec))
    209      end DO
    210      
     245     !DO ll = 1, nlayer + 1
     246     !   zh(:,ll) = H0 * log(pr / (ph(:,ll) + psec))
     247     !end DO
     248           ZH(:,1) = 0.
     249     DO LL = 2, nlayer + 1
     250       !ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
     251        H0bis(:, LL-1) = r * TT(:, LL-1) / g
     252          ZH(:, LL) = ZH(:, LL-1) &
     253           + H0bis(:, LL-1)*(PH(:, LL-1)-PH(:,LL))/PP(:, LL-1)
     254     end DO
     255           ZH(:, 1) = H0 * LOG(PR / (PH(:, 1) + PSEC))
     256
    211257     ! Winds and Brunt Vaisala frequency
    212258     DO ll = 2, nlayer
    … …  
    214260        vh(:, ll) = 0.5 * (vv(:, ll) + vv(:, ll - 1)) ! Meridional wind
    215261
    216         bv(:, ll) = 0.5 * (tt(:, ll) + tt(:, ll - 1)) &
    217                     * r**2 / cpp / H0**2              &
    218                     + (tt(:, ll) - tt(:, ll - 1))     &
    219                     / (zh(:, ll) - zh(:, ll - 1))     &
    220                     * r / H0
    221         bv(:, ll) = sqrt(max(bvsec, bv(:,ll)))
     262       BV(:, LL)= G/(0.5 * (TT(:, LL) + TT(:, LL - 1)))                    &
     263        *((TT(:, LL) - TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1))+        &
     264        G / cpnew(:,LL))
     265
     266       BV(:,LL) =MAX(1.E-12,BV(:,LL)) ! to ensure that BV is positive
     267       BV(:,LL) = MAX(BVSEC,SQRT(BV(:,LL))) ! make sure it is not too small
    222268     end DO
    223269
    … …  
    229275     vh(:, nlayer + 1) = vv(:, nlayer)
    230276
     277!-----------------------------------------------------------------------------------------------------------------------
    231278! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY
    232 !-----------------------------------------
    233 ! The mod function of here a weird arguments
    234 ! are used to produce the waves characteristics
    235 ! in a stochastic way
     279!-----------------------------------------------------------------------------------------------------------------------
     280! The mod function of here a weird arguments are used to produce the waves characteristics in a stochastic way
    236281     DO jw = 1, nw
    237         DO ii = 1, ngrid
     282              !Angle
     283              DO ii = 1, ngrid
    238284           
    239            ran_num_1 = mod(tt(ii, jw) * 10., 1.)
    240            ran_num_2 = mod(tt(ii, jw) * 100., 1.)
    241 
    242            ! angle (0 or pi so far)
    243            zp(jw, ii) = (sign(1., 0.5 - ran_num_1) &
    244                         + 1.) * pi / 2.
    245            ! horizontal wavenumber amplitude
    246            ! TN+GG April/June 2020
    247            ! "Individual waves are not supposed
    248            ! to occupy the entire domain, but
    249            ! only a faction of it" Lott+2012
    250 
    251            !zk(jw, ii) = kmin + (kmax - kmin) * ran_num_2
    252            kstar = pi / sqrt(cell_area(ii))
    253            zk(jw, ii) = max(kmin,kstar) + (kmax - max(kmin,kstar)) * ran_num_2
    254            ! horizontal phase speed
    255            cpha = 0.
    256            DO jj = 1, na
    257               ran_num_3 = mod(tt(ii, jw + 3 * jj)**2, 1.)
    258               cpha = cpha + 2. * cmax *            &
    259                      (ran_num_3 - 0.5) *           &
    260                      sqrt(3.) / sqrt(na * 1.)
    261            end DO
    262            IF (cpha < 0.) THEN
    263               cpha = - 1. * cpha
    264               zp (jw, ii) = zp(jw, ii) + pi
    265            ENDIF
    266            ! Intrinsic frequency
    267            zo(jw, ii) = cpha * zk(jw, ii)
    268            ! Intrinsic frequency is imposed
    269            zo(jw, ii) = zo(jw, ii)                                    &
    270                       + zk(jw, ii) * cos(zp(jw, ii)) * uh(ii, launch) &
     285                ran_num_1 = mod(tt(ii, jw) * 10., 1.)
     286                ran_num_2 = mod(tt(ii, jw) * 100., 1.)
     287
     288                ! Angle (0 or pi so far)
     289                zp(jw, ii) = (sign(1., 0.5 - ran_num_1) + 1.) * pi / 2.
     290                ! Horizontal wavenumber amplitude
     291                ! TN+GG April/June 2020 (rev2381)
     292                ! "Individual waves are not supposed to occupy the entire domain, but only a faction of it" Lott+2012
     293                !zk(jw, ii) = kmin + (kmax - kmin) * ran_num_2
     294                kstar = pi / sqrt(cell_area(ii))
     295                zk(jw, ii) = max(kmin,kstar) + (kmax - max(kmin,kstar)) * ran_num_2
     296               
     297               ! Horizontal phase speed
     298! this computation allows to get a gaussian distribution centered on 0 (right ?)
     299! then cmin is not useful, and it favors small values.
     300                cpha = 0.
     301                DO jj = 1, na
     302                        ran_num_3 = mod(tt(ii, jw + 3 * jj)**2, 1.)
     303                        cpha = cpha + 2. * cmax *                    &
     304                        (ran_num_3 - 0.5) *                          &
     305                        sqrt(3.) / sqrt(na * 1.)
     306                end DO
     307                IF (cpha < 0.) THEN
     308                cpha = - 1. * cpha
     309                zp (jw, ii) = zp(jw, ii) + pi
     310                ENDIF
     311! otherwise, with the computation below, we get a uniform distribution between cmin and cmax.
     312!           ran_num_3 = mod(tt(ii, jw)**2, 1.)
     313!           cpha = cmin + (cmax - cmin) * ran_num_3
     314
     315                ! Intrinsic frequency
     316                zo(jw, ii) = cpha * zk(jw, ii)
     317                ! Intrinsic frequency is imposed
     318                zo(jw, ii) = zo(jw, ii)                                    &
     319                      + zk(jw, ii) * cos(zp(jw, ii)) * uh(ii, launch)      &
    271320                      + zk(jw, ii) * sin(zp(jw, ii)) * vh(ii, launch)
    272            ! Momentum flux at launch level
    273            epflux_0(jw, ii) = epflux_max                              &
     321                ! Momentum flux at launch level
     322                epflux_0(jw, ii) = epflux_max                              &
    274323                        * mod(100. * (uu(ii, jw)**2 + vv(ii, jw)**2), 1.)
    275 
    276         end DO
    277      end DO
     324              ENDDO
     325     end DO !DO jw = 1, nw
    278326
    279327!    print*,'epflux_max just after waves charac. ramdon:', epflux_max
    280328!    print*,'epflux_0 just after waves charac. ramdon:', epflux_0
     329
     330!-----------------------------------------------------------------------------------------------------------------------
    281331! 4. COMPUTE THE FLUXES
    282 !----------------------
    283      ! 4.1 Vertical velocity at launching altitude to ensure
    284      !     the correct value to the imposed fluxes.
     332!-----------------------------------------------------------------------------------------------------------------------
     333     ! 4.1 Vertical velocity at launching altitude to ensure the correct value to the imposed fluxes.
    285334     !------------------------------------------------------
    286335     DO jw = 1, nw
    … …  
    362411                      ! No breaking (Eq. 6):
    363412                      wwm(jw, :)                                                          &
    364                       ! Dissipation (Eq. 8):
    365                       * exp(-rdiss * pr / (ph(:, ll + 1) + ph (:, ll))                    &
     413                      ! Dissipation (Eq. 8)! (real rho used here rather than pressure
     414                      ! parametration because the original code has a bug if the density of
     415                      ! the planet at the launch altitude not approximate 1):             
     416                      * exp(-rdiss * 2./rho(:, LL + 1)                                    &
    366417                      * ((bv(:, ll + 1) + bv (:, ll)) / 2.)**3                            &
    367418                      / max(abs(intr_freq_p(jw, :) + intr_freq_m(jw, :)) / 2., zoisec)**4 &
    … …  
    370421                      ! changes sign)
    371422                      max(0., sign(1., intr_freq_p(jw, :) * intr_freq_m(jw, :)))          &
    372                       ! Saturation (Eq. 12)
    373                       * abs(intr_freq_p(jw, :))**3 / bv(:, ll + 1)                        &
    374                       * exp(-zh(:, ll + 1) / H0)                                          &
     423                      ! Saturation (Eq. 12)(rho at launch altitude is imposed by J.Liu.
     424                      ! Same reason with Eq. 8)
     425                      * abs(intr_freq_p(jw, :))**3 /(2.0* bv(:, ll + 1))                  &
     426                      * rho(:,launch)*exp(-zh(:, ll + 1) / H0)                            &
    375427                      * sat**2 * kmin**2 / zk(jw, :)**4)
    376428        end DO
     429       ! END OF W(KB)ARNING
    377430
    378431     ! Evaluate EP-flux from Eq. 7 and give the right orientation to the stress
    … …  
    388441           v_epflux_tot(:, ll + 1) = v_epflux_tot(:, ll + 1) + v_epflux_p(jw, :)
    389442           east_gwstress(:, ll) = east_gwstress(:, ll)                         & 
    390                                 + max(0., u_epflux_p(jw, :)) / float(nw)
     443                                + max(0., u_epflux_p(jw, ll)) / float(nw)
    391444           west_gwstress(:, ll) = west_gwstress(:, ll)                         &
    392445                                + min(0., u_epflux_p(jw, ll)) / float(nw)
    393446        end DO
    394447     end DO ! DO LL = LAUNCH, nlayer - 1
    395 
    396 !    print*,'u_epflux_tot just after launching:', u_epflux_tot
    397 !    print*,'v_epflux_tot just after launching:', v_epflux_tot
    398 !    print*,'u_epflux_p just after launching:', maxval(u_epflux_p), minval(u_epflux_p)
    399 !    print*,'v_epflux_p just after launching:', maxval(v_epflux_p), minval(v_epflux_p)
    400448
    401449! 5. TENDENCY CALCULATIONS
    … …  
    436484     !---------------------------------------------
    437485     DO LL = 1, nlayer
    438        d_u(:, LL) = (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL))! &
    439        !d_u(:, LL) = G * (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL))! &
    440         !          / (PH(:, LL + 1) - PH(:, LL)) * DTIME
    441        d_v(:, LL) = (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL))! &
    442        !d_v(:, LL) = G * (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL))! &
    443          !         / (PH(:, LL + 1) - PH(:, LL)) * DTIME
     486       !Notice here the d_u and d_v are tendency (i.e. -1/rho*dEP/dz For Mars) but not increment(Venus).
     487       d_u(:, LL) = - (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL)) &
     488                 / (rho(:,ll) * (ZH(:, LL + 1) - ZH(:, LL)))
     489       d_v(:, LL) = - (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL)) &
     490                 / (rho(:,ll) * (ZH(:, LL + 1) - ZH(:, LL)))
    444491     ENDDO
    445492     d_t(:,:) = 0.
    … …  
    454501     dv_nonoro_gwd(:,:) = d_v(:,:)
    455502
    456 !    print*,'u_epflux_tot just after tendency:', u_epflux_tot
    457 !    print*,'v_epflux_tot just after tendency:', v_epflux_tot
    458 !    print*,'d_u just after tendency:', maxval(d_u(:,:)), minval(d_u)
    459 !    print*,'d_v just after tendency:', maxval(d_v(:,:)), minval(d_v)
    460503     ! Cosmetic: evaluation of the cumulated stress
    461 
    462504     ZUSTR(:) = 0.
    463505     ZVSTR(:) = 0.
    464506     DO LL = 1, nlayer
    465         ZUSTR(:) = ZUSTR(:) + D_U(:, LL) !/ g * (PH(:, LL + 1) - PH(:, LL))
    466         ZVSTR(:) = ZVSTR(:) + D_V(:, LL) !/ g * (PH(:, LL + 1) - PH(:, LL))
     507        ZUSTR(:) = ZUSTR(:) + D_U(:, LL) / g * (PH(:, LL + 1) - PH(:, LL))
     508        ZVSTR(:) = ZVSTR(:) + D_V(:, LL) / g * (PH(:, LL + 1) - PH(:, LL))
    467509     ENDDO
    468510

  • trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90

    r2542 r2595  
    6060                              waterrain, global1d, szangle
    6161      use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
    62       use conc_mod
     62      use conc_mod, only: rnew, cpnew, ini_conc_mod
    6363      use phys_state_var_mod
    6464      use callcorrk_mod, only: callcorrk
    … …  
    12691269      IF (calllott_nonoro) THEN
    12701270
    1271          CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep,pplay,  &
     1271         CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep,        &
     1272                    cpnew, rnew,                            &
     1273                    pplay,                                  &
    12721274                    zmax_th,                                &! max altitude reached by thermals (m)
    12731275                    pt, pu, pv,                             &
    … …  
    12791281         pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)      &
    12801282                               + d_t_hin(1:ngrid,1:nlayer)
    1281 !        d_t_hin(:,:)= d_t_hin(:,:)/ptimestep ! K/s (for outputs?)
    12821283         pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer)      &
    12831284                               + d_u_hin(1:ngrid,1:nlayer)
    1284 !        d_u_hin(:,:)= d_u_hin(:,:)/ptimestep ! (m/s)/s (for outputs?)
    12851285         pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer)      &
    12861286                               + d_v_hin(1:ngrid,1:nlayer)
    1287 !        d_v_hin(:,:)= d_v_hin(:,:)/ptimestep ! (m/s)/s (for outputs?)
    12881287         print*,'du_nonoro: max ', maxval(d_u_hin), 'min ', minval(d_u_hin)
    12891288         print*,'dv_nonoro: max ', maxval(d_v_hin), 'min ', minval(d_v_hin)
    … …  
    21762175!        GW non-oro outputs
    21772176      if (calllott_nonoro) then
    2178          call WRITEDIAGFI(ngrid,"dugwno","GW non-oro dU","m/s2", 3,d_u_hin/ptimestep)
    2179          call WRITEDIAGFI(ngrid,"dvgwno","GW non-oro dV","m/s2", 3,d_v_hin/ptimestep)
     2177         call WRITEDIAGFI(ngrid,"dugwno","GW non-oro dU","m/s2", 3,d_u_hin)
     2178         call WRITEDIAGFI(ngrid,"dvgwno","GW non-oro dV","m/s2", 3,d_v_hin)
    21802179      endif
    21812180     
Note: See TracChangeset for help on using the changeset viewer.