Changeset 124 for trunk/libf/dyn3dpar/disvert_noterre.F

Timestamp:

May 19, 2011, 5:05:39 PM (14 years ago)

Author:

emillour

Message:

EM: suite mise au point discretisation verticale et quelques corrections de bugs dans le version de reference parallele.

File:

• trunk/libf/dyn3dpar/disvert_noterre.F (modified) (15 diffs)

trunk/libf/dyn3dpar/disvert_noterre.F

@@ -4 +4 @@
 c    Nouvelle version 100% Mars !!
 c    On l'utilise aussi pour Venus et Titan, legerment modifiee.
+
+#ifdef CPP_IOIPSL
+      use IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+      use ioipsl_getincom
+#endif
 
       IMPLICIT NONE
@@ -12 +19 @@
 #include "comconst.h"
 #include "logic.h"
+#include "iniprint.h"
 c
 c=======================================================================
@@ -24 +32 @@
       INTEGER l,ll
       REAL snorm
-      REAL alpha,beta,gama,delta,deltaz,h,quoi,quand
+      REAL alpha,beta,gama,delta,deltaz
+      real quoi,quand
       REAL zsig(llm),sig(llm+1)
       INTEGER np,ierr
@@ -44 +53 @@
 c-----------------------------------------------------------------------
 c
+! Initializations:
       pi=2.*ASIN(1.)
+      
+      hybrid=.true. ! default value for hybrid (ie: use hybrid coordinates)
+      CALL getin('hybrid',hybrid)
+      write(lunout,*)'disvert_noterre: hybrid=',hybrid
 
 ! Ouverture possible de fichiers typiquement E.T.
@@ -67 +81 @@
 c        <-> energie cinetique, d'apres la note de Frederic Hourdin...
 
-         PRINT*,'*****************************'
-         PRINT*,'WARNING Lecture de esasig.def'
-         PRINT*,'*****************************'
-         READ(99,*) h
+         write(lunout,*)'*****************************'
+         write(lunout,*)'WARNING reading esasig.def'
+         write(lunout,*)'*****************************'
+         READ(99,*) scaleheight
          READ(99,*) dz0
          READ(99,*) dz1
@@ -76 +90 @@
          CLOSE(99)
 
-         dz0=dz0/h
-         dz1=dz1/h
+         dz0=dz0/scaleheight
+         dz1=dz1/scaleheight
 
          sig1=(1.-dz1)/tanh(.5*(llm-1)/nhaut)
@@ -116 +130 @@
 
       ELSE IF(ierr4.eq.0) then
-         PRINT*,'****************************'
-         PRINT*,'Lecture de z2sig.def'
-         PRINT*,'****************************'
-
-         READ(99,*) h
+         write(lunout,*)'****************************'
+         write(lunout,*)'Reading z2sig.def'
+         write(lunout,*)'****************************'
+
+         READ(99,*) scaleheight
          do l=1,llm
             read(99,*) zsig(l)
@@ -128 +142 @@
          sig(1) =1
          do l=2,llm
-           sig(l) = 0.5 * ( exp(-zsig(l)/h) + exp(-zsig(l-1)/h) )
+           sig(l) = 0.5 * ( exp(-zsig(l)/scaleheight) + 
+     &                      exp(-zsig(l-1)/scaleheight) )
          end do
          sig(llm+1) =0
@@ -134 +149 @@
 c-----------------------------------------------------------------------
       ELSE
-         write(*,*) 'didn t you forget something ??? '
-         write(*,*) 'We need the file  z2sig.def ! (OR esasig.def) '
+         write(lunout,*) 'didn t you forget something ??? '
+         write(lunout,*) 'We need file  z2sig.def ! (OR esasig.def)'
          stop
       ENDIF
@@ -156 +171 @@
 c-----------------------------------------------------------------------
 c
-c se trouvait dans Titan... Verifier...
-c     h = 40.
-
-      if (1.eq.1) then  ! toujours hybrides
-         write(*,*) "*******************************"
-         write(*,*) "Systeme en coordonnees hybrides"
-         write(*,*) 
+
+      if (hybrid) then  ! use hybrid coordinates
+         write(lunout,*) "*********************************"
+         write(lunout,*) "Using hybrid vertical coordinates"
+         write(lunout,*) 
 c        Coordonnees hybrides avec mod
          DO l = 1, llm
@@ -172 +185 @@
          bp(llmp1) = 0.
          ap(llmp1) = 0.
-      else
-         write(*,*) "****************************"
-         write(*,*) "Systeme en coordonnees sigma"
-         write(*,*) 
+      else ! use sigma coordinates
+         write(lunout,*) "********************************"
+         write(lunout,*) "Using sigma vertical coordinates"
+         write(lunout,*) 
 c        Pour ne pas passer en coordonnees hybrides
          DO l = 1, llm
@@ -186 +199 @@
       bp(llmp1) =   0.
 
-      PRINT *,' BP '
-      PRINT *,  bp
-      PRINT *,' AP '
-      PRINT *,  ap
+      write(lunout,*)' BP '
+      write(lunout,*)  bp
+      write(lunout,*)' AP '
+      write(lunout,*)  ap
 
 c     Calcul au milieu des couches :
@@ -197 +210 @@
 c     (on met la meme distance (en log pression)  entre P(llm)
 c     et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est
-c     Specifique.  Ce choix est spécifié ici ET dans exner_hyb.F
+c     Specifique.  Ce choix est spécifié ici ET dans exner_milieu.F
 
       DO l = 1, llm-1
@@ -204 +217 @@
       ENDDO
      
-c     if (hybrid) then
+      if (hybrid) then
          aps(llm) = aps(llm-1)**2 / aps(llm-2) 
          bps(llm) = 0.5*(bp(llm) + bp(llm+1))
-c     else
-c        bps(llm) = bps(llm-1)**2 / bps(llm-2) 
-c        aps(llm) = 0.
-c     end if
-
-      PRINT *,' BPs '
-      PRINT *,  bps
-      PRINT *,' APs'
-      PRINT *,  aps
+      else
+         bps(llm) = bps(llm-1)**2 / bps(llm-2) 
+         aps(llm) = 0.
+      end if
+
+      write(lunout,*)' BPs '
+      write(lunout,*)  bps
+      write(lunout,*)' APs'
+      write(lunout,*)  aps
 
       DO l = 1, llm
        presnivs(l) = aps(l)+bps(l)*preff
-       pseudoalt(l) = -h*log(presnivs(l)/preff)
+       pseudoalt(l) = -scaleheight*log(presnivs(l)/preff)
       ENDDO
 
-      PRINT *,' PRESNIVS' 
-      PRINT *,presnivs 
-      PRINT *,'Pseudo altitude des Presnivs : (avec H = ',h,' km)' 
-      PRINT *,pseudoalt
+      write(lunout,*)' PRESNIVS' 
+      write(lunout,*)presnivs 
+      write(lunout,*)'Pseudo altitude of Presnivs : (for a scale ',
+     &                'height of ',scaleheight,' km)' 
+      write(lunout,*)pseudoalt
 
 c     --------------------------------------------------
@@ -232 +246 @@
 c     --------------------------------------------------
 c     open (53,file='testhybrid.tab')
-c     h=15.5
+c     scaleheight=15.5
 c     do iz=0,34
 c       z = -5 + min(iz,34-iz)
 c     approximation of scale height for Venus
-c       h = 15.5 - z/55.*10.
-c       ps = preff*exp(-z/h)
-c       zsig(1)= -h*log((aps(1) + bps(1)*ps)/preff)
+c       scaleheight = 15.5 - z/55.*10.
+c       ps = preff*exp(-z/scaleheight)
+c       zsig(1)= -scaleheight*log((aps(1) + bps(1)*ps)/preff)
 c       do l=2,llm
 c     approximation of scale height for Venus
 c          if (zsig(l-1).le.55.) then
-c             h = 15.5 - zsig(l-1)/55.*10.
+c             scaleheight = 15.5 - zsig(l-1)/55.*10.
 c          else
-c             h = 5.5 - (zsig(l-1)-55.)/35.*2.
+c             scaleheight = 5.5 - (zsig(l-1)-55.)/35.*2.
 c          endif
-c          zsig(l)= zsig(l-1)-h*
+c          zsig(l)= zsig(l-1)-scaleheight*
 c    .    log((aps(l) + bps(l)*ps)/(aps(l-1) + bps(l-1)*ps))
 c       end do