Changeset 127 for trunk/libf/dyn3dpar/disvert.F90

Timestamp:

May 24, 2011, 1:26:29 PM (14 years ago)

Author:

emillour

Message:

Ehouarn: suite de l'implementation de la discretisation verticale,
avec quelques mises a jour pour concorder avec la version terrestre.
-> Finalement, on met un flag "disvert_type" pour fixer la discretisation

disvert_type==1 (defaut si planet_type=="earth") pour cas terrestre
disvert_type==2 (defaut si planet_type!="earth") pour cas planeto (z2sig.def)

-> au passage, pour rester en phase avec modele terrestre on renomme

disvert_terre en disvert (le disvert "alternatif" demeure 'disvert_noterre')

File:

• trunk/libf/dyn3dpar/disvert.F90 (moved) (moved from trunk/libf/dyn3dpar/disvert_terre.F90) (8 diffs)

trunk/libf/dyn3dpar/disvert.F90

@@ -1 @@
-! $Id: disvert.F90 1480 2011-01-31 21:29:58Z jghattas $
+! $Id: disvert.F90 1520 2011-05-23 11:37:09Z emillour $
 
-SUBROUTINE disvert_terre(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, *)'disvert: 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,*)'disvert: ATTENTION discretisation z a ajuster'
+           write(lunout,*) trim(modname), &
+           ' ATTENTION discretisation z a ajuster'
            dsigmin=1.
         endif
-        WRITE(LUNOUT,*) 'disvert: Discretisation verticale DSIGMIN=',dsigmin
+        write(lunout,*) trim(modname), &
+        ' Discretisation verticale DSIGMIN=',dsigmin
      endif
 
@@ -119 +127 @@
   ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
 
-  write(lunout, *)'disvert: BP '
+  write(lunout, *)  trim(modname),': BP '
   write(lunout, *) bp
-  write(lunout, *)'disvert: 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, *) 'disvert: PRESNIVS '
+  write(lunout, *) trim(modname),': PRESNIVS '
   write(lunout, *) presnivs
 
-END SUBROUTINE disvert_terre
+END SUBROUTINE disvert