Changeset 1520 for LMDZ5/trunk/libf/dyn3dpar/disvert.F90

Timestamp:

May 23, 2011, 1:37:09 PM (14 years ago)

Author:

Ehouarn Millour

Message:

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

File:

• LMDZ5/trunk/libf/dyn3dpar/disvert.F90 (modified) (8 diffs)

LMDZ5/trunk/libf/dyn3dpar/disvert.F90

@@ -1 +1 @@
 ! $Id$
 
-SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
+SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight)
 
   ! Auteur : P. Le Van
@@ -17 +17 @@
   ! ds(l) : distance entre les couches l et l-1 en coord.s
 
-  REAL pa, preff
-  REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1)
-  REAL presnivs(llm)
+  real,intent(in) :: pa, preff
+  real,intent(out) :: ap(llmp1), bp(llmp1)
+  real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
+  real,intent(out) :: presnivs(llm)
+  real,intent(out) :: scaleheight
 
   REAL sig(llm+1), dsig(llm)
@@ -26 +28 @@
   INTEGER l
   REAL dsigmin
-  REAL alpha, beta, deltaz, h
+  REAL alpha, beta, deltaz,scaleheight
   INTEGER iostat
-  REAL pi, x
+  REAL x
+  character(len=*),parameter :: modname="disvert"
 
   !-----------------------------------------------------------------------
 
-  pi = 2 * ASIN(1.)
+  ! default scaleheight is 8km for earth
+  scaleheight=8.
 
   OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
@@ -38 +42 @@
   IF (iostat == 0) THEN
      ! cas 1 on lit les options dans sigma.def:
-     READ(99, *) h ! hauteur d'echelle 8.
+     READ(99, *) scaleheight ! hauteur d'echelle 8.
      READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
      READ(99, *) beta ! facteur d'acroissement en haut 1.3
@@ -44 +48 @@
      READ(99, *) k1 ! nombre de couches dans la transition haute
      CLOSE(99)
-     alpha=deltaz/(llm*h)
-     write(lunout, *)'h, alpha, k0, k1, beta'
+     alpha=deltaz/(llm*scaleheight)
+     write(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', &
+                               scaleheight, alpha, k0, k1, beta
 
      alpha=deltaz/tanh(1./k0)*2.
@@ -51 +56 @@
      sig(1)=1.
      do l=1, llm
-        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) &
-             *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
-        zk=-h*log(sig(l+1))
+        sig(l+1)=(cosh(l/k0))**(-alpha*k0/scaleheight) &
+             *exp(-alpha/scaleheight*tanh((llm-k1)/k0) &
+                  *beta**(l-(llm-k1))/log(beta))
+        zk=-scaleheight*log(sig(l+1))
 
         dzk1=alpha*tanh(l/k0)
@@ -73 +79 @@
            dsigmin=1.
         else
-           WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
+           write(lunout,*) trim(modname), &
+           ' ATTENTION discretisation z a ajuster'
            dsigmin=1.
         endif
-        WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
+        write(lunout,*) trim(modname), &
+        ' Discretisation verticale DSIGMIN=',dsigmin
      endif
 
@@ -119 +127 @@
   ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
 
-  write(lunout, *)' BP '
+  write(lunout, *)  trim(modname),': BP '
   write(lunout, *) bp
-  write(lunout, *)' AP '
+  write(lunout, *)  trim(modname),': AP '
   write(lunout, *) ap
 
   write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
   write(lunout, *)'couches calcules pour une pression de surface =', preff
-  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de'
-  write(lunout, *)'8km'
+  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de '
+  write(lunout, *) scaleheight,' km'
   DO l = 1, llm
      dpres(l) = bp(l) - bp(l+1)
      presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
      write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
-          log(preff/presnivs(l))*8. &
-          , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ &
+          log(preff/presnivs(l))*scaleheight &
+          , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
           max(ap(l+1)+bp(l+1)*preff, 1.e-10))
   ENDDO
 
-  write(lunout, *) 'PRESNIVS '
+  write(lunout, *) trim(modname),': PRESNIVS '
   write(lunout, *) presnivs