Ignore:
Timestamp:
Feb 11, 2019, 5:53:10 PM (6 years ago)
Author:
jvatant
Message:

+ Fix some ambiguities about effective gravity field. Choice is made to keep coherence between dynamics and physics by default,

but having correct altitudes for chemistry - to have correct fluxes - no matter what's in physics.
Then chemistry is "recast" on the dynphy grid.

+ Also correct the fact that chemistry requires altitutde with ref to the geoid not to the local surface -> pphi+phhis in the calculation.
+ Both points above lead to new zlay_eff/zlev_eff altitudes variables in physiq_mod sent to chemistry.

We still need to think about it for microphysics.

--JVO

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

Legend:

Unmodified
Added
Removed

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

    r2053 r2097  
    11SUBROUTINE calchim(ngrid,qy_c,declin,dtchim,            &
    2      ctemp,cpphi,cplay,cplev,czlay,czlev,dqyc)
     2     ctemp,cpphi,cpphis,cplay,cplev,czlay,czlev,dqyc)
    33
    44  !---------------------------------------------------------------------------------------------------------
    … …  
    3131  !                          ( and there'll be work to do to get rid of averaged fields )
    3232  !
     33  !               + 02/19 : To always have correct photodissociations rates, altitudes sent here by physiq are always
     34  !                         calculated with effective g - and with reference to the body not the local surface -
     35  !                          even if in physiq we keep altitudes coherent with dynamics !
     36  !
    3337  ! + STILL TO DO : Replug the interaction with haze (cf titan.old) -> to see with JB.
    3438  !---------------------------------------------------------------------------------------------------------
    … …  
    8084  REAL*8, DIMENSION(ngrid,klev),      INTENT(IN)   :: ctemp       ! Mid-layer temperature (K).
    8185  REAL*8, DIMENSION(ngrid,klev),      INTENT(IN)   :: cpphi       ! Mid-layer geopotential (m2.s-2).
     86  REAL*8, DIMENSION(ngrid),           INTENT(IN)   :: cpphis      ! Surface geopotential (m2.s-2).
    8287  REAL*8, DIMENSION(ngrid,klev),      INTENT(IN)   :: cplay       ! Mid-layer pressure (Pa).
    8388  REAL*8, DIMENSION(ngrid,klev+1),    INTENT(IN)   :: cplev       ! Inter-layer pressure (Pa).
    84   REAL*8, DIMENSION(ngrid,klev),      INTENT(IN)   :: czlay       ! Mid-layer altitude (m).
    85   REAL*8, DIMENSION(ngrid,klev+1),    INTENT(IN)   :: czlev       ! Inter-layer altitude (m).
     89  REAL*8, DIMENSION(ngrid,klev),      INTENT(IN)   :: czlay       ! Mid-layer effective altitude (m) : ref = geoid.
     90  REAL*8, DIMENSION(ngrid,klev+1),    INTENT(IN)   :: czlev       ! Inter-layer effective altitude (m) ref = geoid.
    8691
    8792  REAL*8, DIMENSION(ngrid,klev,nkim), INTENT(OUT)  :: dqyc        ! Chemical species tendencies on GCM layers (mol/mol/s).
    … …  
    201206     ! JVO18 : altitudes are no more calculated in firstcall, as I set kedd in pressure grid
    202207
    203      ! a. For GCM layers we just copy-paste
     208     ! a. For GCM layers we just copy-paste ( assuming that physiq always send correct altitudes ! )
    204209
    205210     PRINT*,'Init chemistry : pressure, density, temperature ... :'
    … …  
    223228
    224229     ! b. Extension in upper atmosphere with Vervack profile
    225      !    NB : Maybe the transition klev/klev+1 is harsh ...     
     230     ! NB : Maybe the transition klev/klev+1 is harsh if T profile different from Vervack ...     
    226231
    227232     ipres=1
    … …  
    270275     !      Below GCM top we have dynamic mixing and for levs < nld=klev-15 the chem. solver ignores diffusion
    271276
    272      ! First calculate kedd for upper chemistry layers
     277     !! First calculate kedd for upper chemistry layers
     278     !DO l=klev-4,nlaykim_tot
     279     !   logp=-log10(press_c(l))
     280     !! 2E6 at 400 km ~ 10-2 mbar
     281     !   IF     ( logp.ge.2.0 .and. logp.le.3.0 ) THEN
     282     !         kedd(l) = 2.e6 * 5.0**(logp-2.0)
     283     !! 1E7 at 500 km ~ 10-3 mbar
     284     !   ELSE IF     ( logp.ge.3.0 .and. logp.le.4.0 ) THEN
     285     !         kedd(l) = 1.e7 * 3.0**(logp-3.0)
     286     !! 3E7 above 700 km ~ 10-4 mbar
     287     !   ELSEIF ( logp.gt.4.0                   ) THEN
     288     !        kedd(l) = 3.e7
     289     !   ENDIF
     290     !ENDDO
     291
     292     ! Kedd from (E7) in Vuitton 2019
    273293     DO l=klev-4,nlaykim_tot
    274         logp=-log10(press_c(l))
    275      ! 2E6 at 400 km ~ 10-2 mbar
    276         IF     ( logp.ge.2.0 .and. logp.le.3.0 ) THEN
    277               kedd(l) = 2.e6 * 5.0**(logp-2.0)
    278      ! 1E7 at 500 km ~ 10-3 mbar
    279         ELSE IF     ( logp.ge.3.0 .and. logp.le.4.0 ) THEN
    280               kedd(l) = 1.e7 * 3.0**(logp-3.0)
    281      ! 3E7 above 700 km ~ 10-4 mbar
    282         ELSEIF ( logp.gt.4.0                   ) THEN
    283              kedd(l) = 3.e7
    284         ENDIF
     294       kedd(l) = 300.0 * ( 1.0E2 / press_c(l) )**1.5 * 3.0E7 /  &
     295               ( 300.0 * ( 1.0E2 / press_c(l) )**1.5 + 3.0E7 )
    285296     ENDDO
    286297     
    … …  
    323334
    324335        ! a. For GCM layers we just copy-paste
     336        ! JVO 19 : Now physiq always sent correct altitudes with effective g for chemistry ( even if it's not the case in physiq )
    325337
    326338        DO l=1,klev
    … …  
    335347
    336348        ! b. Extension in upper atmosphere with Vervack profile
    337         !    NB : The transition klev/klev+1 is maybe harsh if we have correct g above but not under 
    338349
    339350        ipres=1
    … …  
    356367        DO l=klev+1,nlaykim_tot
    357368
    358            ! Compute geopotential on the upper grid, to have correct altitude
    359            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    360            ! NB : Upper chemistry model is based on observational profiles and needs correct altitudes !!
    361            ! Altitudes can (and must) be calculated with correct g even if it's not the case in the GCM under
     369           ! Compute geopotential on the upper grid with effective g to have correct altitudes
     370           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    362371
    363372           temp1 = 0.5*(temp_c(l-1)+temp_c(l)) ! interlayer temp
    364373           phi_c(l) = phi_c(l-1) + r*temp1*log(press_c(l-1)/press_c(l)) ! Geopotential assuming hydrostatic equilibrium
    365374
    366            rmil(l)  =  ( g*rad*rad / (g*rad - phi_c(l)) ) / 1000.0 ! z(phi) with g varying with altitude
     375           rmil(l)  =  ( g*rad*rad / (g*rad - ( phi_c(l) + cpphis(ig) ) ) ) / 1000.0 ! z(phi) with g varying with altitude with reference to the geoid
    367376        ENDDO
    368377
    369378        DO l=klev+2,nlaykim_tot
    370           rinter(l) = 0.5*(rmil(l-1) + rmil(l))
     379          rinter(l) = 0.5*(rmil(l-1) + rmil(l)) ! should be balanced with the intermediate pressure rather than 0.5
    371380        ENDDO
    372381

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

    r2050 r2097  
    260260      real zlss                      ! Sub solar point longitude (radians).
    261261      real zday                      ! Date (time since Ls=0, calculated in sols).
    262       real zzlay(ngrid,nlayer)       ! Altitude at the middle of the atmospheric layers.
    263       real zzlev(ngrid,nlayer+1)     ! Altitude at the atmospheric layer boundaries.
    264 
     262      real zzlay(ngrid,nlayer)       ! Altitude at the middle of the atmospheric layers (ref : local surf).
     263      real zzlev(ngrid,nlayer+1)     ! Altitude at the atmospheric layer boundaries (ref : local surf).
     264      real zzlay_eff(ngrid,nlayer)   ! Effective altitude at the middle of the atmospheric layers (ref : geoid ).
     265      real zzlev_eff(ngrid,nlayer+1) ! Effective altitude at the atmospheric layer boundaries ( ref : geoid ).
     266     
    265267! TENDENCIES due to various processes :
    266268
    … …  
    704706     
    705707        do l=1,nlayer         
    706            zzlay(:,l) = pphi(:,l) / gzlat(:,l)
     708           zzlay(:,l) = pphi(:,l) / gzlat(:,l) ! Reference = local surface
    707709        enddo
    708710       
    … …  
    710712     
    711713        do l=1,nlayer
    712            zzlay(:,l) = g*rad*rad / ( g*rad - pphi(:,l) ) - rad
     714           zzlay(:,l) = g*rad*rad / ( g*rad - ( pphi(:,l) + phisfi(:) )) - rad
    713715           gzlat(:,l) = g*rad*rad / ( rad + zzlay(:,l) )**2
    714716        end do
    … …  
    718720      zzlev(:,1)=0.
    719721      zzlev(:,nlayer+1)=1.e7 ! Dummy top of last layer above 10000 km...
     722      ! JVO 19 : This altitude is indeed dummy for the GCM and fits ptop=0
     723      ! but for upper chemistry that's a pb -> we anyway redefine it just after ..
    720724
    721725      do l=2,nlayer
    … …  
    727731      enddo     
    728732
    729       ! Zonal averages needed for chemistry
    730       ! ~~~~~~~ Taken from old Titan ~~~~~~
    731       if (moyzon_ch) then
    732          
    733          if (eff_gz .eqv. .false.) then
    734            zzlaybar(1,:)=(zphibar(1,:)+zphisbar(1))/g
    735          else ! if we take into account effective altitude-dependent gravity field
    736            zzlaybar(1,:)=g*rad*rad/(g*rad-(zphibar(1,:)+zphisbar(1)))-rad
    737          endif
     733      ! Effective altitudes ( eg needed for chemistry ) with correct g, and with reference to the geoid
     734      ! JVO 19 : We shall always have correct altitudes in chemistry no matter what's in physics
     735      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     736      if (moyzon_ch) then ! Zonal averages
     737         
     738         zzlaybar(1,:)=g*rad*rad/(g*rad-(zphibar(1,:)+zphisbar(1)))-rad   ! reference = geoid
    738739         zzlevbar(1,1)=zphisbar(1)/g       
    739740         DO l=2,nlayer
    … …  
    747748            if (latitude(ig).ne.latitude(ig-1)) then       
    748749               DO l=1,nlayer   
    749                   if (eff_gz .eqv. .false.) then
    750                     zzlaybar(ig,l)=(zphibar(ig,l)+zphisbar(ig))/g
    751                   else ! if we take into account effective altitude-dependent gravity field
    752                     zzlaybar(ig,l)=g*rad*rad/(g*rad-(zphibar(ig,l)+zphisbar(ig)))-rad
    753                   endif           
     750                  zzlaybar(ig,l)=g*rad*rad/(g*rad-(zphibar(ig,l)+zphisbar(ig)))-rad
    754751               ENDDO
    755752               zzlevbar(ig,1)=zphisbar(ig)/g
    … …  
    765762            endif         
    766763         ENDDO
     764
     765      else !  if not moyzon
     766     
     767         DO l=1,nlayer   
     768           zzlay_eff(ig,l)=g*rad*rad/(g*rad-(pphi(ig,l)+phisfi(ig)))-rad ! reference = geoid
     769         ENDDO
     770         zzlev_eff(ig,1)=phisfi(ig)/g
     771         DO l=2,nlayer
     772            z1=(pplay(ig,l-1)+pplev(ig,l))/ (pplay(ig,l-1)-pplev(ig,l))
     773            z2=(pplev(ig,l)  +pplay(ig,l))/(pplev(ig,l)  -pplay(ig,l))
     774            zzlev_eff(ig,l)=(z1*zzlay_eff(ig,l-1)+z2*zzlay_eff(ig,l))/(z1+z2)
     775         ENDDO
     776         zzlev_eff(ig,nlayer+1)=zzlay_eff(ig,nlayer)+(zzlay_eff(ig,nlayer)-zzlev_eff(ig,nlayer))
    767777
    768778      endif  ! moyzon
    … …  
    10751085            tpq(:,:,:) = tpq(:,:,:) + zdqmufi(:,:,:)*ptimestep ! only manipulation of tpq->*ptimestep here
    10761086#else
    1077             call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,zdqmufi)
     1087            call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,zdqmufi) ! JVO 19 : To be fixed, what altitude do we need ?
    10781088            pdq(:,:,:) = pdq(:,:,:) + zdqmufi(:,:,:)
    10791089#endif
    … …  
    11731183
    11741184                 ! Here we send zonal average fields ( corrected with cond ) from dynamics to chem. module
    1175                  call calchim(ngrid,ychimbar,declin,ctimestep,ztfibar,zphibar,  &
     1185                 call calchim(ngrid,ychimbar,declin,ctimestep,ztfibar,zphibar,zphisbar,  &
    11761186                              zplaybar,zplevbar,zzlaybar,zzlevbar,dycchi)
    11771187               else ! 3D chemistry (or 1D run)
    1178                  call calchim(ngrid,ychim,declin,ctimestep,pt,pphi,  &
    1179                               pplay,pplev,zzlay,zzlev,dycchi)
     1188                 call calchim(ngrid,ychim,declin,ctimestep,pt,pphi,phisfi,  &
     1189                              pplay,pplev,zzlay_eff,zzlev_eff,dycchi)
    11801190               endif ! if moyzon
    11811191
    … …  
    16021612      CALL send_xios_field("dtrad",dtrad)
    16031613      CALL send_xios_field("dtdyn",zdtdyn)
     1614      CALL send_xios_field("dtdif",zdtdif)
    16041615
    16051616      ! Chemical tracers
Note: See TracChangeset for help on using the changeset viewer.