Changeset 1194 for trunk/LMDZ.GENERIC/libf

Timestamp:

Mar 4, 2014, 11:32:02 AM (11 years ago)

Author:

sglmd

Message:

Latitude-dependent gravity field added. Option oblate = .true. in callphys.def, and three additional variables needed: J2, Rmean and MassPlanet?.

Location:

trunk/LMDZ.GENERIC/libf/phystd

Files:

• callcorrk.F90 (modified) (4 diffs)
• callkeys.h (modified) (4 diffs)
• inifis.F (modified) (1 diff)
• mass_redistribution.F90 (modified) (2 diffs)
• optci.F90 (modified) (2 diffs)
• optcv.F90 (modified) (2 diffs)
• physiq.F90 (modified) (5 diffs)
• radcommon_h.F90 (modified) (1 diff)
• turbdiff.F90 (modified) (2 diffs)

trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90

@@ -180 +180 @@
       REAL*8 muvarrad(L_LEVELS)
 
-         
 !===============================================================
 !     Initialization on first call
@@ -227 +226 @@
          call suaer_corrk       ! set up aerosol optical properties
 
-         Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed
 
          if((igcm_h2o_vap.eq.0) .and. varactive)then
@@ -714 +712 @@
 
 
+         Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed
+         glat_ig=glat(ig)
+
 !-----------------------------------------------------------------------
 !     Shortwave
 
          if(fract(ig) .ge. 1.0e-4) then ! only during daylight!
-!            if(sum(eclipse) .eq. 0.) then
-!                   write(*,*) 'but eclipse is still 0 in physiq !'
-!               endif
             if((ngrid.eq.1).and.(global1d))then
                do nw=1,L_NSPECTV
@@ -801 +799 @@
          DO l=2,L_NLAYRAD
             dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))  &
-                *g/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
+                *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
             dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))  &
-                *g/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
+                *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
          END DO      
 
 !     These are values at top of atmosphere
          dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv)           &
-             *g/(cpp*scalep*(plevrad(3)-plevrad(1)))
+             *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(1)))
          dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi)           &
-             *g/(cpp*scalep*(plevrad(3)-plevrad(1)))
+             *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(1)))
 
 !     ---------------------------------------------------------------

trunk/LMDZ.GENERIC/libf/phystd/callkeys.h

@@ -12 +12 @@
      &   , meanOLR, specOLR                                             &
      &   , kastprof, Tstrat                                             &
-     &   , nosurf, intheat, flatten                                     &     
+     &   , nosurf, intheat, flatten, oblate                             &     
      &   , newtonian, tau_relax, testradtimes                           &
      &   , check_cpp_match, force_cpp                                   &
@@ -39 +39 @@
      &   , cloudlvl                                                     &
      &   , pceil                                                        &
+     &   , Rmean, J2, MassPlanet                                        &
      &   , strictboundcorrk                                                                                     
 
@@ -74 +75 @@
       logical CLFvarying
       logical nosurf
+      logical oblate
 
       integer iddist
@@ -108 +110 @@
       real intheat
       real flatten
+      real Rmean
+      real J2
+      real MassPlanet

trunk/LMDZ.GENERIC/libf/phystd/inifis.F

@@ -158 +158 @@
          write(*,*) "rings_shadow = ", rings_shadow
          
-         write(*,*) "Flattening of the planet (a-b)/a ?"
+         write(*,*) "Compute latitude-dependent gravity field?"
+         oblate = .false.
+         call getin("oblate", oblate)
+         write(*,*) "oblate = ", oblate
+
+         write(*,*) "Flattening of the planet (a-b)/a "
          flatten = 0.0
          call getin("flatten", flatten)
          write(*,*) "flatten = ", flatten
+         
+
+         write(*,*) "Needed if oblate=.true.: J2"
+         J2 = 0.0
+         call getin("J2", J2)
+         write(*,*) "J2 = ", J2
+         
+         write(*,*) "Needed if oblate=.true.: Planet mass (*1e24 kg)"
+         MassPlanet = 0.0
+         call getin("MassPlanet", MassPlanet)
+         write(*,*) "MassPlanet = ", MassPlanet         
+
+         write(*,*) "Needed if oblate=.true.: Planet mean radius (m)"
+         Rmean = 0.0
+         call getin("Rmean", Rmean)
+         write(*,*) "Rmean = ", Rmean
          
 ! Test of incompatibility:

trunk/LMDZ.GENERIC/libf/phystd/mass_redistribution.F90

@@ -6 +6 @@
        USE watercommon_h, Only: Tsat_water,RLVTT
        use surfdat_h
+       use radcommon_h, only: glat
        USE comgeomfi_h
        USE tracer_h
@@ -128 +129 @@
          zdmass_sum(ig,nlayermx+1)=0.
          DO l = nlayermx, 1, -1
-           zmass(ig,l) = (pplev(ig,l)-pplev(ig,l+1))/g
+           zmass(ig,l) = (pplev(ig,l)-pplev(ig,l+1))/glat(ig)
            zdmass_sum(ig,l)= zdmass_sum(ig,l+1)+pdmassmr(ig,l)
          END DO

trunk/LMDZ.GENERIC/libf/phystd/optci.F90

@@ -4 +4 @@
 
   use radinc_h
-  use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi
+  use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig
   use gases_h
   implicit none
@@ -123 +123 @@
         U(k)  = Cmk*DPR(k)*mugaz/muvar(k) 
      else
-        dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))*mugaz/muvar(k)
+        dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k)
         U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in optci.F  
             !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present

trunk/LMDZ.GENERIC/libf/phystd/optcv.F90

@@ -4 +4 @@
 
   use radinc_h
-  use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep,indv
+  use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep,indv,glat_ig
   use gases_h
 
@@ -111 +111 @@
         U(k)  = Cmk*DPR(k)*mugaz/muvar(k) 
      else
-        dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))*mugaz/muvar(k)
+        dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k)
         U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in optci.F  
             !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present

trunk/LMDZ.GENERIC/libf/phystd/physiq.F90

@@ -11 +11 @@
       use watercommon_h, only : RLVTT, Psat_water,epsi
       use gases_h
-      use radcommon_h, only: sigma, eclipse
+      use radcommon_h, only: sigma, eclipse, glat, grav
       use radii_mod, only: h2o_reffrad, co2_reffrad, h2o_cloudrad
       use aerosol_mod
@@ -267 +267 @@
       real ztime_fin
       real zdh(ngrid,nlayermx)
+      real gmplanet
+
       integer length
       parameter (length=100)
@@ -471 +473 @@
          !! this is defined in radcommon_h
          ALLOCATE(eclipse(ngrid))
+         ALLOCATE(glat(ngrid))
 
 !        variables set to 0
@@ -648 +651 @@
       end if
 
+
+
+!    Compute variations of g with latitude (oblate case)
+!
+        if (oblate .eq. .false.) then
+           glat(:) = g
+        else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then
+        print*,'I need values for flatten, J2, Rmean and MassPlanet to compute glat (else set oblate=.false.)'
+        call abort
+        else
+
+        gmplanet = MassPlanet*grav*1e24
+        do ig=1,ngrid
+           glat(ig)= gmplanet/(Rmean**2) * (1.D0 + 0.75 *J2 - 2.0*flatten/3. + (2.*flatten - 15./4.* J2) * cos(2. * (pi/2. - lati(ig)))) 
+        end do
+        endif
+
+!!      write(*,*) 'lati, glat =',lati(1)*180./pi,glat(1), lati(ngrid/2)*180./pi, glat(ngrid/2)
+
+
+
 !     Compute geopotential between layers
 !     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-      zzlay(1:ngrid,1:nlayermx)=pphi(1:ngrid,1:nlayermx)/g
+      zzlay(1:ngrid,1:nlayermx)=pphi(1:ngrid,1:nlayermx)
+      do l=1,nlayermx         
+      zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
+      enddo
 
       zzlev(1:ngrid,1)=0.
@@ -672 +699 @@
             zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
             zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
-            mass(ig,l)  = (pplev(ig,l) - pplev(ig,l+1))/g
+            mass(ig,l)  = (pplev(ig,l) - pplev(ig,l+1))/glat(ig)
             massarea(ig,l)=mass(ig,l)*area(ig)
          enddo
       enddo
+
 
 !-----------------------------------------------------------------------

trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90

@@ -131 +131 @@
 
       real*8, parameter :: sigma = 5.67032D-8
+      real*8, parameter :: grav = 6.672E-11
 
       real*8 Cmk
       save Cmk
+      real*8 glat_ig
+      save glat_ig
 
       ! extinction of incoming sunlight (Saturn's rings, eclipses, etc...)
       REAL, DIMENSION(:), ALLOCATABLE :: eclipse
 
+      !Latitude-dependent gravity
+      REAL, DIMENSION(:), ALLOCATABLE :: glat
+
       end module radcommon_h

trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90

@@ -8 +8 @@
 
       use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol,Psat_water
-      use radcommon_h, only : sigma
+      use radcommon_h, only : sigma, glat
       USE surfdat_h
       USE comgeomfi_h
@@ -165 +165 @@
       DO ilay=1,nlay
          DO ig=1,ngrid
-            zmass(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
+            zmass(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/glat(ig)
             zExner(ig,ilay)=(pplev(ig,ilay)/pplev(ig,1))**rcp
             zovExner(ig,ilay)=1./ppopsk(ig,ilay)