Changeset 1947 for trunk/LMDZ.TITAN/libf


Ignore:
Timestamp:
Jun 13, 2018, 5:22:44 PM (7 years ago)
Author:
jvatant
Message:

Enables altitude-depending gravity field g(z) (glat->gzlat) in physics
+ Can be dangerous ( disagreement with dyn) but important (compulsive !) to have correct altitudes in the chemistry
+ Can be activated with eff_gz flag in callphys.def
-- JVO

Location:
trunk/LMDZ.TITAN/libf/phytitan
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.TITAN/libf/phytitan/calchim.F90

    r1946 r1947  
    6262  USE geometry_mod, ONLY: latitude
    6363  USE logic_mod, ONLY: moyzon_ch
    64   USE moyzon_mod, ONLY: tmoy, playmoy, phimoy, zlaymoy, zlevmoy
     64  USE moyzon_mod, ONLY: tmoy, playmoy
    6565
    6666  IMPLICIT NONE
    … …  
    198198     ! ----------------------------------------------------------------------
    199199     
    200      ! JVO18 : altitudes calculated in firstcall are just for diag, as I set kedd in pressure grid
     200     ! JVO18 : altitudes are no more calculated in firstcall, as I set kedd in pressure grid
    201201
    202202     ! a. For GCM layers we just copy-paste
    203203
    204 !$OMP MASTER
    205204     PRINT*,'Init chemistry : pressure, density, temperature ... :'
    206      PRINT*,'level, rmil (km), press_c (mbar), nb (cm-3), temp_c(K)'
    207 !$OMP END MASTER
    208 !$OMP BARRIER
     205     PRINT*,'level, press_c (mbar), nb (cm-3), temp_c (K)'
    209206     
    210207     IF (ngrid.NE.1) THEN
    211208       DO l=1,klev
    212           rinter(l)  = (zlevmoy(l)+rad)/1000.0              ! km
    213           rmil(l)    = (zlaymoy(l)+rad)/1000.0              ! km
    214209          temp_c(l)  = tmoy(l)                              ! K
    215           phi_c(l)   = phimoy(l)                            ! m2.s-2
    216210          press_c(l) = playmoy(l)/100.                      ! mbar
    217211          nb(l)      = 1.e-4*press_c(l) / (kbol*temp_c(l))  ! cm-3
    218 !$OMP MASTER
    219           PRINT*, l, rmil(l)-rad/1000.0, press_c(l), nb(l), temp_c(l), phi_c(l)
    220 !$OMP END MASTER
    221 !$OMP BARRIER
     212          PRINT*, l, press_c(l), nb(l), temp_c(l), phi_c(l)
    222213       ENDDO
    223        rinter(klev+1)=(zlevmoy(klev+1)+rad)/1000.
    224214     ELSE
    225215       DO l=1,klev
    226           rinter(l)  = (czlev(1,l)+rad)/1000.0                ! km
    227           rmil(l)    = (czlay(1,l)+rad)/1000.0                ! km
    228216          temp_c(l)  = ctemp(1,l)                             ! K
    229           phi_c(l)   = cpphi(1,l)                             ! m2.s-2
    230217          press_c(l) = cplay(1,l)/100.                        ! mbar
    231218          nb(l)      = 1.e-4*press_c(l) / (kbol*temp_c(l))  ! cm-3
    232           PRINT*, l, rmil(l)-rad/1000.0, press_c(l), nb(l), temp_c(l), phi_c(l)
     219          PRINT*, l, press_c(l), nb(l), temp_c(l)
    233220       ENDDO
    234        rinter(klev+1)=(czlev(1,klev+1)+rad)/1000.
    235      ENDIF ! if ngrid.ne.1
     221     ENDIF
    236222
    237223     ! b. Extension in upper atmosphere with Vervack profile
    … …  
    250236        nb(l)      = exp( log(ct1d(ipres))*(1-factp) + log(ct1d(ipres+1))* factp )
    251237        temp_c(l)  = t1d(ipres)*(1-factp) + t1d(ipres+1)*factp
    252      ENDDO
    253 
    254      ! We build altitude with hydrostatic equilibrium on preskim grid with Vervack profile
    255      ! ( keeping in mind that preskim is built based on Vervack profile with dz=10km )
    256 
    257      DO l=klev+1,nlaykim_tot
    258 
    259         ! Compute geopotential on the upper grid, to have correct altitude
    260         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    261         ! NB : Upper chemistry model is based on observational profiles and needs correct altitudes !!
    262         ! Altitudes can (and must) be calculated with correct g even if it's not the case in the GCM under
    263 
    264         temp1 = 0.5*(temp_c(l-1)+temp_c(l)) ! interlayer temp
    265         phi_c(l) = phi_c(l-1) + r*temp1*log(press_c(l-1)/press_c(l)) ! Geopotential assuming hydrostatic equilibrium
    266      
    267         rmil(l)  =  ( g*rad*rad / (g*rad - phi_c(l)) ) / 1000.0 ! z(phi) with g varying with altitude
    268 
    269 !$OMP MASTER
    270         PRINT*, l , rmil(l)-rad/1000.0, press_c(l), nb(l), temp_c(l), phi_c(l)
    271 !$OMP END MASTER
    272 !$OMP BARRIER
    273      ENDDO
    274 
    275      DO l=klev+2,nlaykim_tot
    276         rinter(l) = 0.5*(rmil(l-1) + rmil(l))
     238        PRINT*, l , press_c(l), nb(l), temp_c(l)
    277239     ENDDO
    278240

  • trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90

    r1897 r1947  
    279279!=======================================================================
    280280
    281 
    282          Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
    283          glat_ig=glat(ig)
     281         Cmk(:)      = 0.01 * 1.0 / (gzlat(ig,:) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
     282         gzlat_ig(:) = gzlat(ig,:)
    284283
    285284!-----------------------------------------------------------------------
    … …  
    386385         DO l=2,L_NLAYRAD
    387386            dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))  &
    388                 *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
     387                *gzlat(ig,L_NLAYRAD+1-l)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
    389388            dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))  &
    390                 *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
     389                *gzlat(ig,L_NLAYRAD+1-l)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
    391390         END DO     
    392391
    393392!     These are values at top of atmosphere
    394393         dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv)           &
    395              *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(1)))
     394             *gzlat(ig,L_NLAYRAD)/(cpp*scalep*(plevrad(3)-plevrad(1)))
    396395         dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi)           &
    397              *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(1)))
     396             *gzlat(ig,L_NLAYRAD)/(cpp*scalep*(plevrad(3)-plevrad(1)))
    398397
    399398

  • trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90

    r1897 r1947  
    4040      logical,save :: oblate
    4141!$OMP THREADPRIVATE(nosurf,oblate)
     42      logical,save :: eff_gz
     43!$OMP THREADPRIVATE(eff_gz)
    4244     
    4345      integer,save :: ichim

  • trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90

    r1926 r1947  
    11
    22
    3 SUBROUTINE calmufi(dt, plev, zlev, play, zlay, temp, pq, zdqfi, zdq)
     3SUBROUTINE calmufi(dt, plev, zlev, play, zlay, g3d, temp, pq, zdqfi, zdq)
    44  !! Interface subroutine to YAMMS model for Titan LMDZ GCM.
    55  !!
    … …  
    1818  USE MMP_GCM
    1919  USE tracer_h
    20   USE comcstfi_mod, only : g
    2120  USE callkeys_mod, only : callclouds
    2221  USE muphy_diag
    … …  
    2928  REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: play  !! Pressure layers (Pa).
    3029  REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlay  !! Altitude at the center of each layer (m).
     30  REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: g3d   !! Latitude-Altitude depending gravitational acceleration (m.s-2).
    3131  REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp  !! Temperature at the center of each layer (K).
    3232
    … …  
    5656  REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: dgazs !! Tendencies of each condensible gaz species !(\(mol.mol^{-1}\)).
    5757
    58   REAL(kind=8), DIMENSION(:), ALLOCATABLE ::  int2ext
     58  REAL(kind=8), DIMENSION(:,:), ALLOCATABLE ::  int2ext
    5959  TYPE(error) :: err
    6060
    … …  
    6868  nices = size(ices_indx)
    6969  ! Conversion intensive to extensive
    70   ALLOCATE( int2ext(nlay) ) 
     70  ALLOCATE( int2ext(nlon,nlay) ) 
    7171
    7272  ! Allocate arrays
    … …  
    101101    ! Convert tracers to extensive ( except for gazs where we work with molar mass ratio )
    102102    ! We suppose a given order of tracers !
    103     int2ext(:) = ( plev(ilon,1:nlay)-plev(ilon,2:nlay+1) ) / g
    104     m0as(:) = zq(ilon,:,1) * int2ext(:)
    105     m3as(:) = zq(ilon,:,2) * int2ext(:)
    106     m0af(:) = zq(ilon,:,3) * int2ext(:)
    107     m3af(:) = zq(ilon,:,4) * int2ext(:)
     103    int2ext(ilon,:) = ( plev(ilon,1:nlay)-plev(ilon,2:nlay+1) ) / g3d(ilon,1:nlay)
     104    m0as(:) = zq(ilon,:,1) * int2ext(ilon,:)
     105    m3as(:) = zq(ilon,:,2) * int2ext(ilon,:)
     106    m0af(:) = zq(ilon,:,3) * int2ext(ilon,:)
     107    m3af(:) = zq(ilon,:,4) * int2ext(ilon,:)
    108108   
    109109    if (callclouds) then ! if call clouds
    110       m0n(:) = zq(ilon,:,5) * int2ext(:)
    111       m3n(:) = zq(ilon,:,6) * int2ext(:)
     110      m0n(:) = zq(ilon,:,5) * int2ext(ilon,:)
     111      m3n(:) = zq(ilon,:,6) * int2ext(ilon,:)
    112112      do i=1,nices
    113         m3i(:,nices) = zq(ilon,:,6+i) * int2ext(:)
     113        m3i(:,nices) = zq(ilon,:,6+i) * int2ext(ilon,:)
    114114        gazs(:,i)    = zq(ilon,:,ices_indx(i)) / rat_mmol(ices_indx(i)) ! For gazs we work on the full tracer array !!
    115115        ! We use the molar mass ratio from GCM in case there is discrepancy with the mm one
    … …  
    151151    ! We suppose a given order of tracers !
    152152 
    153     zdq(ilon,:,1) = dm0as(:) / int2ext(:)
    154     zdq(ilon,:,2) = dm3as(:) / int2ext(:)
    155     zdq(ilon,:,3) = dm0af(:) / int2ext(:)
    156     zdq(ilon,:,4) = dm3af(:) / int2ext(:)
     153    zdq(ilon,:,1) = dm0as(:) / int2ext(ilon,:)
     154    zdq(ilon,:,2) = dm3as(:) / int2ext(ilon,:)
     155    zdq(ilon,:,3) = dm0af(:) / int2ext(ilon,:)
     156    zdq(ilon,:,4) = dm3af(:) / int2ext(ilon,:)
    157157   
    158158    if (callclouds) then ! if call clouds
    159       zdq(ilon,:,5) = dm0n(:) / int2ext(:)
    160       zdq(ilon,:,6) = dm3n(:) / int2ext(:)
     159      zdq(ilon,:,5) = dm0n(:) / int2ext(ilon,:)
     160      zdq(ilon,:,6) = dm3n(:) / int2ext(ilon,:)
    161161      do i=1,nices
    162         zdq(ilon,:,6+i) = dm3i(:,nices) / int2ext(:)
     162        zdq(ilon,:,6+i) = dm3i(:,nices) / int2ext(ilon,:)
    163163        zdq(ilon,:,ices_indx(i)) = dgazs(:,i) * rat_mmol(ices_indx(i)) ! For gazs we work on the full tracer array !!
    164164        ! We use the molar mass ratio from GCM in case there is discrepancy with the mm one

  • trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90

    r1897 r1947  
    181181     call getin_p("Rmean", Rmean)
    182182     write(*,*) "Rmean = ", Rmean
     183     
     184     write(*,*) "Compute effective altitude-dependent gravity field?"
     185     eff_gz = .false.
     186     call getin_p("eff_gz", eff_gz)
     187     write(*,*) "eff_gz = ", eff_gz
    183188         
    184189! Test of incompatibility:
    … …  
    344349       stop
    345350     endif
     351     
     352     ! sanity check warning
     353     if (callchim.and.(.not.eff_gz)) then
     354       print*,"WARNING : You run chemistry without effective altitude-dependent gravity field !!"
     355       print*,"You will have wrong vertical photodissociations rates, that's crazy !!"
     356       print*,"I let you continue but you should rather set eff_gz =.true. ..."
     357       !stop
     358     endif
     359
    346360
    347361     write(*,*)"Chemistry is computed every ichim", &

  • trunk/LMDZ.TITAN/libf/phytitan/inimufi.F90

    r1926 r1947  
    22
    33  use mmp_gcm
    4   use callkeys_mod, only : callclouds, p_prod, tx_prod, rc_prod, air_rad
     4  use callkeys_mod, only : callclouds, p_prod, tx_prod, rc_prod, air_rad, eff_gz
    55  use tracer_h
    66  use comcstfi_mod, only : g, rad, mugaz
    … …  
    4040
    4141  ! PATCH : YAMMS now allows to enable/disable effective g computations:
    42   !  In the library it defaults to .true., in the GCM we do not want it !
    43   mm_use_effg = .false.
     42  mm_use_effg = eff_gz
    4443
    4544  !----------------------------------------------------

  • trunk/LMDZ.TITAN/libf/phytitan/mass_redistribution.F90

    r1647 r1947  
    55                                                   
    66       use surfdat_h
    7        use radcommon_h, only: glat
     7       use radcommon_h, only: gzlat
    88       USE tracer_h
    99       USE planete_mod, only: bp
    … …  
    123123         zdmass_sum(ig,nlayer+1)=0.
    124124         DO l = nlayer, 1, -1
    125            zmass(ig,l) = (pplev(ig,l)-pplev(ig,l+1))/glat(ig)
     125           zmass(ig,l) = (pplev(ig,l)-pplev(ig,l+1))/gzlat(ig,l)
    126126           zdmass_sum(ig,l)= zdmass_sum(ig,l+1)+pdmassmr(ig,l)
    127127         END DO

  • trunk/LMDZ.TITAN/libf/phytitan/optci.F90

    r1897 r1947  
    33
    44  use radinc_h
    5   use radcommon_h, only: gasi,tlimit,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig,gweight
     5  use radcommon_h, only: gasi,tlimit,Cmk,gzlat_ig,tgasref,pfgasref,wnoi,scalep,indi,gweight
    66  use gases_h
    7   use comcstfi_mod, only: g, r
     7  use comcstfi_mod, only: r
    88  use callkeys_mod, only: continuum,graybody,callclouds,callmufi, uncoupl_optic_haze
    99  use tracer_h, only : nmicro,nice
    … …  
    118118
    119119  do K=2,L_LEVELS
     120 
     121     ilay = k / 2 ! int. arithmetic => gives the gcm layer index
     122     
    120123     DPR(k) = PLEV(K)-PLEV(K-1)
    121124
    122125     ! if we have continuum opacities, we need dz
    123       dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))
    124       U(k)  = Cmk*DPR(k)     ! only Cmk line in optci.F
     126      dz(k) = dpr(k)*R*TMID(K)/(gzlat_ig(ilay)*PMID(K))
     127      U(k)  = Cmk(ilay)*DPR(k)     ! only Cmk line in optci.F
    125128     
    126129     call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))
    … …  
    221224           enddo
    222225
    223            ! Oobleck test
    224            !rho = PMID(k)*scalep / (TMID(k)*286.99)
    225            !if(WNOI(nw).gt.300.0 .and. WNOI(nw).lt.500.0)then
    226            !   DCONT = rho * 0.125 * 4.6e-4
    227            !elseif(WNOI(nw).gt.500.0 .and. WNOI(nw).lt.700.0)then
    228            !   DCONT = 1000*dpr(k) * 1.0 * 4.6e-4 / g
    229            !   DCONT = rho * 1.0 * 4.6e-4
    230            !elseif(WNOI(nw).gt.700.0 .and. WNOI(nw).lt.900.0)then
    231            !   DCONT = rho * 0.125 * 4.6e-4
    232            !endif
    233 
    234226           DCONT = DCONT*dz(k)
    235227

  • trunk/LMDZ.TITAN/libf/phytitan/optcv.F90

    r1897 r1947  
    33
    44  use radinc_h
    5   use radcommon_h, only: gasv, tlimit, Cmk, tgasref, pfgasref,wnov,scalep,indv,glat_ig,gweight
     5  use radcommon_h, only: gasv,tlimit,Cmk,gzlat_ig,tgasref,pfgasref,wnov,scalep,indv,gweight
    66  use gases_h
    7   use comcstfi_mod, only: g, r
     7  use comcstfi_mod, only: r
    88  use callkeys_mod, only: continuum,graybody,callgasvis,callclouds,callmufi,uncoupl_optic_haze
    99  use tracer_h, only: nmicro,nice
    … …  
    132132
    133133  do K=2,L_LEVELS
     134 
     135     ilay = k / 2 ! int. arithmetic => gives the gcm layer index
     136 
    134137     DPR(k) = PLEV(K)-PLEV(K-1)
    135138
    136139     ! if we have continuum opacities, we need dz
    137140
    138       dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))
    139       U(k)  = Cmk*DPR(k)     ! only Cmk line in optcv.F     
     141      dz(k) = dpr(k)*R*TMID(K)/(gzlat_ig(ilay)*PMID(K))
     142      U(k)  = Cmk(ilay)*DPR(k)     ! only Cmk line in optcv.F     
    140143
    141144     call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))

  • trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90

    r1943 r1947  
    1515 
    1616      use radinc_h, only : L_NSPECTI,L_NSPECTV
    17       use radcommon_h, only: sigma, glat, grav, BWNV
     17      use radcommon_h, only: sigma, gzlat, gzlat_ig, Cmk, grav, BWNV
    1818      use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe
    1919      use comchem_h, only: nkim, cnames, nlaykim_up, ykim_up, ykim_tot
    … …  
    394394      integer :: i,j
    395395      real :: pqo(ngrid,nlayer,nq)   ! Tracers for the optics (X/m2).
    396       real :: i2e(nlayer)            ! int 2 ext factor
    397396      ! -----******----- END FOR MUPHYS OPTICS -----******-----
    398397
    399       real,dimension(:,:), allocatable :: i2e2d ! factor to convert X.kg-1 in X.m-3 (for microphysics diags)
    400      
     398      real :: i2e(ngrid,nlayer)      ! int 2 ext factor ( X.kg-1 -> X.m-3 or X.m-2 , depending if optics or diags )
     399
    401400      real,save,dimension(:,:,:), allocatable :: tpq ! Tracers for decoupled microphysical tests ( temporary in 01/18 )
    402401!$OMP THREADPRIVATE(tpq)
    … …  
    410409    !   Or one can put calmufi in MMP_GCM module (in muphytitan).
    411410    INTERFACE
    412       SUBROUTINE calmufi(dt, plev, zlev, play, zlay, temp, pq, zdqfi, zdq)
     411      SUBROUTINE calmufi(dt, plev, zlev, play, zlay, g3d, temp, pq, zdqfi, zdq)
    413412        REAL(kind=8), INTENT(IN)                 :: dt    !! Physics timestep (s).
    414413        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: plev  !! Pressure levels (Pa).
    … …  
    416415        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: play  !! Pressure layers (Pa).
    417416        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlay  !! Altitude at the center of each layer (m).
     417        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: g3d   !! Latitude-Altitude depending gravitational acceleration (m.s-2).
    418418        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp  !! Temperature at the center of each layer (K). 
    419419        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: pq    !! Tracers (\(kg.kg^{-1}}\)).
    … …  
    476476        ALLOCATE(fract(ngrid))           
    477477         ! This is defined in radcommon_h
    478         ALLOCATE(glat(ngrid))
    479        
     478        ALLOCATE(gzlat(ngrid,nlayer))
     479        ALLOCATE(gzlat_ig(nlayer))
     480        ALLOCATE(Cmk(nlayer))         
    480481
    481482!        Variables set to 0
    … …  
    704705      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    705706      if (oblate .eqv. .false.) then     
    706          glat(:) = g         
     707         gzlat(:,:) = g         
    707708      else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then     
    708          print*,'I need values for flatten, J2, Rmean and MassPlanet to compute glat (else set oblate=.false.)'
     709         print*,'I need values for flatten, J2, Rmean and MassPlanet to compute gzlat (else set oblate=.false.)'
    709710         call abort
    710711      else
    711712         gmplanet = MassPlanet*grav*1e24
    712713         do ig=1,ngrid
    713             glat(ig)= gmplanet/(Rmean**2) * (1.D0 + 0.75 *J2 - 2.0*flatten/3. + (2.*flatten - 15./4.* J2) * cos(2. * (pi/2. - latitude(ig))))
     714            gzlat(ig,:)= gmplanet/(Rmean**2) * (1.D0 + 0.75 *J2 - 2.0*flatten/3. + (2.*flatten - 15./4.* J2) * cos(2. * (pi/2. - latitude(ig))))
    714715         end do
    715716      endif
    716 
    717       ! Compute geopotential between layers.
    718       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    719       zzlay(:,:)=pphi(:,:)
    720       do l=1,nlayer         
    721          zzlay(:,l)= zzlay(:,l)/glat(:)
    722       enddo
     717     
     718      ! Compute altitudes with the geopotential coming from the dynamics.
     719      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     720
     721      if (eff_gz .eqv. .false.) then
     722     
     723        do l=1,nlayer         
     724           zzlay(:,l) = pphi(:,l) / gzlat(:,l)
     725        enddo
     726       
     727      else ! In this case we consider variations of g with altitude
     728     
     729        do l=1,nlayer
     730           zzlay(:,l) = g*rad*rad / ( g*rad - pphi(:,l) ) - rad
     731           gzlat(:,l) = g*rad*rad / ( rad + zzlay(:,l) )**2
     732        end do
     733     
     734      endif ! if eff_gz
    723735
    724736      zzlev(:,1)=0.
    … …  
    737749      if (moyzon_ch) then
    738750         
    739          zzlaybar(1,:)=(zphibar(1,:)+zphisbar(1))/g
    740          ! SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
    741          !       zzlaybar(1,:)=RG*RA*RA/(RG*RA-(zphibar(1,:)+zphisbar(1)))-RA
    742          zzlevbar(1,1)=zphisbar(1)/g
     751         if (eff_gz .eqv. .false.) then
     752           zzlaybar(1,:)=(zphibar(1,:)+zphisbar(1))/g
     753         else ! if we take into account effective altitude-dependent gravity field
     754           zzlaybar(1,:)=g*rad*rad/(g*rad-(zphibar(1,:)+zphisbar(1)))-rad
     755         endif
     756         zzlevbar(1,1)=zphisbar(1)/g       
    743757         DO l=2,nlayer
    744             z1=(zplaybar(1,l-1)+zplevbar(1,l))/(zplevbar(1,l-1)-zplevbar(1,l))
     758            z1=(zplaybar(1,l-1)+zplevbar(1,l))/(zplaybar(1,l-1)-zplevbar(1,l))
    745759            z2=(zplevbar(1,l)  +zplaybar(1,l))/(zplevbar(1,l)  -zplaybar(1,l))
    746760            zzlevbar(1,l)=(z1*zzlaybar(1,l-1)+z2*zzlaybar(1,l))/(z1+z2)
    … …  
    749763
    750764         DO ig=2,ngrid
    751             if (latitude(ig).ne.latitude(ig-1)) then
    752                DO l=1,nlayer
    753                   zzlaybar(ig,l)=(zphibar(ig,l)+zphisbar(ig))/g
    754                   ! SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
    755                   !zzlaybar(ig,l)=RG*RA*RA/(RG*RA-(zphibar(ig,l)+zphisbar(ig)))-RA
     765            if (latitude(ig).ne.latitude(ig-1)) then       
     766               DO l=1,nlayer   
     767                  if (eff_gz .eqv. .false.) then
     768                    zzlaybar(ig,l)=(zphibar(ig,l)+zphisbar(ig))/g
     769                  else ! if we take into account effective altitude-dependent gravity field
     770                    zzlaybar(ig,l)=g*rad*rad/(g*rad-(zphibar(ig,l)+zphisbar(ig)))-rad
     771                  endif           
    756772               ENDDO
    757773               zzlevbar(ig,1)=zphisbar(ig)/g
    758774               DO l=2,nlayer
    759                   z1=(zplaybar(ig,l-1)+zplevbar(ig,l))/ (zplevbar(ig,l-1)-zplevbar(ig,l))
     775                  z1=(zplaybar(ig,l-1)+zplevbar(ig,l))/ (zplaybar(ig,l-1)-zplevbar(ig,l))
    760776                  z2=(zplevbar(ig,l)  +zplaybar(ig,l))/(zplevbar(ig,l)  -zplaybar(ig,l))
    761777                  zzlevbar(ig,l)=(z1*zzlaybar(ig,l-1)+z2*zzlaybar(ig,l))/(z1+z2)
    762778               ENDDO
    763                zzlevbar(ig,nlayer+1)=zzlaybar(ig,nlayer)+(zzlaybar(ig,nlayer)-zzlevbar(ig,nlayer))
     779               zzlevbar(ig,nlayer+1)=zzlaybar(ig,nlayer)+(zzlaybar(ig,nlayer)-zzlevbar(ig,nlayer))         
    764780            else
    765781               zzlaybar(ig,:)=zzlaybar(ig-1,:)
    766782               zzlevbar(ig,:)=zzlevbar(ig-1,:)
    767             endif
     783            endif         
    768784         ENDDO
    769785
    … …  
    778794            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
    779795            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
    780             mass(ig,l)  = (pplev(ig,l) - pplev(ig,l+1))/glat(ig)
     796            mass(ig,l)  = (pplev(ig,l) - pplev(ig,l+1))/gzlat(ig,l)
    781797            massarea(ig,l)=mass(ig,l)*cell_area(ig)
    782798         enddo
    … …  
    843859               ! NOTE: it should be moved somewhere else: calmufi performs the same kind of
    844860               ! computations... waste of time...
    845                DO i = 1, ngrid
    846                  i2e(1:nlayer) = ( pplev(i,1:nlayer)-pplev(i,2:nlayer+1) ) / g
    847                  pqo(i,:,:) = 0.0
    848                  DO j=1,nmicro-nice
    849                    pqo(i,:,j) = pq(i,:,j)*i2e(:)
    850                  ENDDO
     861               i2e(:,1:nlayer) = ( pplev(:,1:nlayer)-pplev(:,2:nlayer+1) ) / gzlat(:,1:nlayer)
     862               pqo(:,:,:) = 0.0
     863               DO j=1,nmicro-nice
     864                 pqo(:,:,j) = pq(:,:,j)*i2e(:,:)
    851865               ENDDO
    852 
    853866
    854867               ! standard callcorrk
    … …  
    10861099#ifdef USE_QTEST
    10871100            dtpq(:,:,:) = 0.0 ! we want tpq to go only through mufi
    1088             call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,pt,tpq,dtpq,zdqmufi)
     1101            call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,tpq,dtpq,zdqmufi)
    10891102            tpq(:,:,:) = tpq(:,:,:) + zdqmufi(:,:,:)*ptimestep ! only manipulation of tpq->*ptimestep here
    10901103#else
    1091             call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,pt,pq,pdq,zdqmufi)
     1104            call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,zdqmufi)
    10921105            pdq(:,:,:) = pdq(:,:,:) + zdqmufi(:,:,:)
    10931106#endif
    … …  
    10971110              zdqfibar(:,:,:) = 0.0 ! We work in zonal average -> forget processes other than condensation
    10981111              call calmufi(ptimestep,zplevbar,zzlevbar,zplaybar,zzlaybar, &
    1099                            ztfibar,zqfibar,zdqfibar,zdqmufibar)
     1112                           gzlat,ztfibar,zqfibar,zdqfibar,zdqmufibar)
    11001113            endif
    11011114
    … …  
    15111524         if (callmufi) then
    15121525           ! Microphysical tracers are expressed in unit/m3.
    1513            ! convert X.kg-1 --> X.m-3
    1514            ALLOCATE(i2e2d(ngrid,nlayer))
    1515            i2e2d(:,:) = ( pplev(:,1:nlayer)-pplev(:,2:nlayer+1) ) / g /(zzlev(:,2:nlayer+1)-zzlev(:,1:nlayer))
     1526           ! convert X.kg-1 --> X.m-3 (whereas for optics was -> X.m-2)
     1527           i2e(:,:) = ( pplev(:,1:nlayer)-pplev(:,2:nlayer+1) ) / gzlat(:,1:nlayer) /(zzlev(:,2:nlayer+1)-zzlev(:,1:nlayer))
    15161528
    15171529#ifdef USE_QTEST
    15181530            ! Microphysical tracers passed through dyn+phys(except mufi)
    1519             call writediagfi(ngrid,"mu_m0as_dp","Dynphys only spherical mode 0th order moment",'m-3',3,zq(:,:,micro_indx(1))*i2e2d)
    1520             call writediagfi(ngrid,"mu_m3as_dp","Dynphys only spherical mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(2))*i2e2d)
    1521             call writediagfi(ngrid,"mu_m0af_dp","Dynphys only fractal mode 0th order moment",'m-3',3,zq(:,:,micro_indx(3))*i2e2d)
    1522             call writediagfi(ngrid,"mu_m3af_dp","Dynphys only fractal mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(4))*i2e2d)
     1531            call writediagfi(ngrid,"mu_m0as_dp","Dynphys only spherical mode 0th order moment",'m-3',3,zq(:,:,micro_indx(1))*i2e)
     1532            call writediagfi(ngrid,"mu_m3as_dp","Dynphys only spherical mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(2))*i2e)
     1533            call writediagfi(ngrid,"mu_m0af_dp","Dynphys only fractal mode 0th order moment",'m-3',3,zq(:,:,micro_indx(3))*i2e)
     1534            call writediagfi(ngrid,"mu_m3af_dp","Dynphys only fractal mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(4))*i2e)
    15231535            ! Microphysical tracers passed through mufi only
    1524             call writediagfi(ngrid,"mu_m0as_mo","Mufi only spherical mode 0th order moment",'m-3',3,tpq(:,:,micro_indx(1))*i2e2d)
    1525             call writediagfi(ngrid,"mu_m3as_mo","Mufi only spherical mode 3rd order moment",'m3/m3',3,tpq(:,:,micro_indx(2))*i2e2d)
    1526             call writediagfi(ngrid,"mu_m0af_mo","Mufi only fractal mode 0th order moment",'m-3',3,tpq(:,:,micro_indx(3))*i2e2d)
    1527             call writediagfi(ngrid,"mu_m3af_mo","Mufi only fractal mode 3rd order moment",'m3/m3',3,tpq(:,:,micro_indx(4))*i2e2d)
     1536            call writediagfi(ngrid,"mu_m0as_mo","Mufi only spherical mode 0th order moment",'m-3',3,tpq(:,:,micro_indx(1))*i2e)
     1537            call writediagfi(ngrid,"mu_m3as_mo","Mufi only spherical mode 3rd order moment",'m3/m3',3,tpq(:,:,micro_indx(2))*i2e)
     1538            call writediagfi(ngrid,"mu_m0af_mo","Mufi only fractal mode 0th order moment",'m-3',3,tpq(:,:,micro_indx(3))*i2e)
     1539            call writediagfi(ngrid,"mu_m3af_mo","Mufi only fractal mode 3rd order moment",'m3/m3',3,tpq(:,:,micro_indx(4))*i2e)
    15281540#else
    1529             call writediagfi(ngrid,"mu_m0as","Spherical mode 0th order moment",'m-3',3,zq(:,:,micro_indx(1))*i2e2d)
    1530             call writediagfi(ngrid,"mu_m3as","Spherical mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(2))*i2e2d)
    1531             call writediagfi(ngrid,"mu_m0af","Fractal mode 0th order moment",'m-3',3,zq(:,:,micro_indx(3))*i2e2d)
    1532             call writediagfi(ngrid,"mu_m3af","Fractal mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(4))*i2e2d)
     1541            call writediagfi(ngrid,"mu_m0as","Spherical mode 0th order moment",'m-3',3,zq(:,:,micro_indx(1))*i2e)
     1542            call writediagfi(ngrid,"mu_m3as","Spherical mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(2))*i2e)
     1543            call writediagfi(ngrid,"mu_m0af","Fractal mode 0th order moment",'m-3',3,zq(:,:,micro_indx(3))*i2e)
     1544            call writediagfi(ngrid,"mu_m3af","Fractal mode 3rd order moment",'m3/m3',3,zq(:,:,micro_indx(4))*i2e)
    15331545#endif
    15341546           
    1535             DEALLOCATE(i2e2d)
    1536 
    15371547            ! Microphysical diagnostics
    15381548            call writediagfi(ngrid,"mmd_aer_prec","Total aerosols precipitations",'m',2,mmd_aer_prec)
    … …  
    16191629         
    16201630        enddo
    1621        
     1631       
     1632        ! Condensation tendencies ( mol/mol/s )
     1633        CALL send_xios_field("dqcond_CH4",zdqcond(:,:,7+nmicro)/rat_mmol(7+nmicro))
     1634        CALL send_xios_field("dqcond_C2H2",zdqcond(:,:,10+nmicro)/rat_mmol(10+nmicro))
     1635        CALL send_xios_field("dqcond_C2H4",zdqcond(:,:,12+nmicro)/rat_mmol(12+nmicro))
     1636        CALL send_xios_field("dqcond_C2H6",zdqcond(:,:,14+nmicro)/rat_mmol(14+nmicro))
     1637        CALL send_xios_field("dqcond_C3H6",zdqcond(:,:,17+nmicro)/rat_mmol(17+nmicro))
     1638        CALL send_xios_field("dqcond_C4H4",zdqcond(:,:,21+nmicro)/rat_mmol(21+nmicro))
     1639        CALL send_xios_field("dqcond_CH3CCH",zdqcond(:,:,24+nmicro)/rat_mmol(24+nmicro))
     1640        CALL send_xios_field("dqcond_C3H8",zdqcond(:,:,25+nmicro)/rat_mmol(25+nmicro))
     1641        CALL send_xios_field("dqcond_C4H2",zdqcond(:,:,26+nmicro)/rat_mmol(26+nmicro))
     1642        CALL send_xios_field("dqcond_C4H6",zdqcond(:,:,27+nmicro)/rat_mmol(27+nmicro))
     1643        CALL send_xios_field("dqcond_C4H10",zdqcond(:,:,28+nmicro)/rat_mmol(28+nmicro))
     1644        CALL send_xios_field("dqcond_AC6H6",zdqcond(:,:,29+nmicro)/rat_mmol(29+nmicro))
     1645        CALL send_xios_field("dqcond_HCN",zdqcond(:,:,36+nmicro)/rat_mmol(36+nmicro))
     1646        CALL send_xios_field("dqcond_CH3CN",zdqcond(:,:,40+nmicro)/rat_mmol(40+nmicro))
     1647        CALL send_xios_field("dqcond_HC3N",zdqcond(:,:,42+nmicro)/rat_mmol(42+nmicro))
     1648        CALL send_xios_field("dqcond_NCCN",zdqcond(:,:,43+nmicro)/rat_mmol(43+nmicro))
     1649        CALL send_xios_field("dqcond_C4N2",zdqcond(:,:,44+nmicro)/rat_mmol(44+nmicro))
     1650
    16221651      endif ! of 'if callchim'
    16231652
    16241653      ! Microphysical tracers
    16251654      if (callmufi) then
    1626         CALL send_xios_field("mu_m0as",zq(:,:,micro_indx(1))*i2e2d)
    1627         CALL send_xios_field("mu_m3as",zq(:,:,micro_indx(2))*i2e2d)
    1628         CALL send_xios_field("mu_m0af",zq(:,:,micro_indx(3))*i2e2d)
    1629         CALL send_xios_field("mu_m3af",zq(:,:,micro_indx(4))*i2e2d)
     1655        CALL send_xios_field("mu_m0as",zq(:,:,micro_indx(1))*i2e)
     1656        CALL send_xios_field("mu_m3as",zq(:,:,micro_indx(2))*i2e)
     1657        CALL send_xios_field("mu_m0af",zq(:,:,micro_indx(3))*i2e)
     1658        CALL send_xios_field("mu_m3af",zq(:,:,micro_indx(4))*i2e)
    16301659      endif       
    16311660

  • trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90

    r1788 r1947  
    8888      real*8, parameter :: grav = 6.672E-11
    8989
    90       real*8,save :: Cmk
    91       real*8,save :: glat_ig
    92 !$OMP THREADPRIVATE(Cmk,glat_ig)
    93 
    9490      ! extinction of incoming sunlight (Saturn's rings, eclipses, etc...)
    9591      REAL, DIMENSION(:), ALLOCATABLE ,SAVE :: eclipse
     92!$OMP THREADPRIVATE(eclipse)
    9693
    97       !Latitude-dependent gravity
    98       REAL, DIMENSION(:), ALLOCATABLE , SAVE :: glat
    99 !$OMP THREADPRIVATE(glat,eclipse)
     94      ! Altitude-Latitude-dependent gravity
     95      REAL, DIMENSION(:,:), ALLOCATABLE , SAVE :: gzlat ! This should be stored elsewhere ...
     96      REAL, DIMENSION(:), ALLOCATABLE , SAVE :: gzlat_ig
     97      REAL, DIMENSION(:), ALLOCATABLE , SAVE :: Cmk
     98!$OMP THREADPRIVATE(gzlat,gzlat_ig,Cmk)
    10099
    101100end module radcommon_h

  • trunk/LMDZ.TITAN/libf/phytitan/turbdiff.F90

    r1647 r1947  
    77          pdqdif,pdqsdif,flux_u,flux_v,lastcall)
    88
    9       use radcommon_h, only : sigma, glat
     9      use radcommon_h, only : sigma, gzlat
    1010      use comcstfi_mod, only: rcp, g, r, cpp
    1111      use callkeys_mod, only: tracer,nosurf
    … …  
    139139      DO ilay=1,nlay
    140140         DO ig=1,ngrid
    141             zmass(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/glat(ig)
     141            zmass(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/gzlat(ig,ilay)
    142142            zExner(ig,ilay)=(pplev(ig,ilay)/pplev(ig,1))**rcp
    143143            zovExner(ig,ilay)=1./ppopsk(ig,ilay)
Note: See TracChangeset for help on using the changeset viewer.