Changeset 2040 for LMDZ5/trunk/libf/dyn3d_common/disvert.F90

Timestamp:

May 9, 2014, 12:50:34 PM (11 years ago)

Author:

Ehouarn Millour

Message:

Added a "strato_custom" mode in disvert (for development purposes).
This is triggered by setting "vert_sampling=strato_custom" in a def file.

More work is needed to clean disvert (and merge it with disvert_noterre).
DC+EM

File:

• LMDZ5/trunk/libf/dyn3d_common/disvert.F90 (modified) (6 diffs)

LMDZ5/trunk/libf/dyn3d_common/disvert.F90

@@ -1 +1 @@
 ! $Id$
 
-SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight)
-
-  ! Auteur : P. Le Van
-
+SUBROUTINE disvert()
+
+#ifdef CPP_IOIPSL
+  use ioipsl, only: getin
+#else
+  USE ioipsl_getincom, only: getin
+#endif
   use new_unit_m, only: new_unit
-  use ioipsl, only: getin
   use assert_m, only: assert
 
@@ -13 +15 @@
   include "dimensions.h"
   include "paramet.h"
+  include "comvert.h"
+  include "comconst.h"
   include "iniprint.h"
   include "logic.h"
 
-  ! s = sigma ** kappa : coordonnee verticale
-  ! dsig(l) : epaisseur de la couche l ds la coord. s
-  ! sig(l) : sigma a l'interface des couches l et l-1
-  ! ds(l) : distance entre les couches l et l-1 en coord.s
-
-  real,intent(in) :: pa, preff
-  real,intent(out) :: ap(llmp1) ! in Pa
-  real, intent(out):: bp(llmp1)
-  real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
-  real,intent(out) :: presnivs(llm)
-  real,intent(out) :: scaleheight
-
+!-------------------------------------------------------------------------------
+! Purpose: Vertical distribution functions for LMDZ.
+!          Triggered by the levels number llm.
+!-------------------------------------------------------------------------------
+! Read    in "comvert.h":
+! pa                         !--- PURE PRESSURE COORDINATE FOR P<pa (in Pascals)
+! preff                      !--- REFERENCE PRESSURE                 (101325 Pa)
+! Written in "comvert.h":
+! ap(llm+1), bp(llm+1)       !--- Ap, Bp HYBRID COEFFICIENTS AT INTERFACES
+! aps(llm),  bps(llm)        !--- Ap, Bp HYBRID COEFFICIENTS AT MID-LAYERS
+! dpres(llm)                 !--- PRESSURE DIFFERENCE FOR EACH LAYER
+! presnivs(llm)              !--- PRESSURE AT EACH MID-LAYER
+! scaleheight                !--- VERTICAL SCALE HEIGHT            (Earth: 8kms)
+! nivsig(llm+1)              !--- SIGMA INDEX OF EACH LAYER INTERFACE
+! nivsigs(llm)               !--- SIGMA INDEX OF EACH MID-LAYER
+!-------------------------------------------------------------------------------
+! Local variables:
   REAL sig(llm+1), dsig(llm)
-  real zk, zkm1, dzk1, dzk2, k0, k1
+  REAL sig0(llm+1), zz(llm+1)
+  REAL zk, zkm1, dzk1, dzk2, z, k0, k1
 
   INTEGER l, unit
@@ -37 +47 @@
   character(len=*),parameter :: modname="disvert"
 
-  character(len=6):: vert_sampling
+  character(len=24):: vert_sampling
   ! (allowed values are "param", "tropo", "strato" and "read")
 
   !-----------------------------------------------------------------------
 
-  print *, "Call sequence information: disvert"
+  WRITE(lunout,*) TRIM(modname)//" starts"
 
   ! default scaleheight is 8km for earth
@@ -49 +59 @@
   vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
   call getin('vert_sampling', vert_sampling)
-  print *, 'vert_sampling = ' // vert_sampling
+  WRITE(lunout,*) TRIM(modname)//' vert_sampling = ' // vert_sampling
   if (llm==39 .and. vert_sampling=="strato") then
      dsigmin=0.3 ! Vieille option par défaut pour CMIP5
@@ -144 +154 @@
      ap(1)=0.
      ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
+  CASE("strato_custom") 
+    ! custumize strato distribution with specific alpha & beta values and function
+    ! depending on llm (experimental and temporary)!
+    SELECT CASE (llm)
+      CASE(55)
+        alpha=0.45
+        beta=4.0
+      CASE(63)
+        alpha=0.45
+        beta=5.0
+      CASE(71)
+        alpha=3.05
+        beta=65.
+      CASE(79)
+        alpha=3.20
+        ! alpha=2.05 ! FLOTT 79 (PLANTE)
+        beta=70.
+    END SELECT
+    ! Or used values provided by user in def file:
+    CALL getin("strato_alpha",alpha)
+    CALL getin("strato_beta",beta)
+    
+    ! Build geometrical distribution
+    scaleheight=7.
+    zz(1)=0.
+    IF (llm==55.OR.llm==63) THEN
+      DO l=1,llm
+        z=zz(l)/scaleheight
+        zz(l+1)=zz(l)+0.03+z*1.5*(1.-TANH(z-0.5))+alpha*(1.+TANH(z-1.5))     &
+                            +5.0*EXP((l-llm)/beta)
+      ENDDO
+    ELSEIF (llm==71) THEN !.OR.llm==79) THEN      ! FLOTT 79 (PLANTE)
+      DO l=1,llm
+        z=zz(l)
+        zz(l+1)=zz(l)+0.02+0.88*TANH(z/2.5)+alpha*(1.+TANH((z-beta)/15.))
+      ENDDO
+    ELSEIF (llm==79) THEN
+      DO l=1,llm
+        z=zz(l)
+        zz(l+1)=zz(l)+0.02+0.80*TANH(z/3.8)+alpha*(1+TANH((z-beta)/17.))     &
+                            +0.03*TANH(z/.25)
+      ENDDO
+    ENDIF ! of IF (llm==55.OR.llm==63) ...
+    
+    ! Build sigma distribution
+    sig0=EXP(-zz(:)/scaleheight)
+    sig0(llm+1)=0.
+    sig=ridders(sig0)
+    
+    ! Compute ap() and bp()
+    bp(1)=1.
+    bp(2:llm)=EXP(1.-1./sig(2:llm)**2)
+    bp(llm+1)=0.
+    ap=pa*(sig-bp)
+
   case("read")
      ! Read "ap" and "bp". First line is skipped (title line). "ap"
@@ -191 +256 @@
   write(lunout, *) presnivs
 
+CONTAINS
+
+!-------------------------------------------------------------------------------
+!
+FUNCTION ridders(sig) RESULT(sg)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Purpose: Search for s solving (Pa/Preff)*s+(1-Pa/Preff)*EXP(1-1./s**2)=sg
+! Notes:   Uses Ridders' method, quite robust. Initial bracketing: 0<=sg<=1.
+! Reference: Ridders, C. F. J. "A New Algorithm for Computing a Single Root of a
+!       Real Continuous Function" IEEE Trans. Circuits Systems 26, 979-980, 1979
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)  :: sig(:)
+  REAL              :: sg(SIZE(sig))
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: it, ns, maxit
+  REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
+!-------------------------------------------------------------------------------
+  ns=SIZE(sig); maxit=9999
+  c1=Pa/Preff; c2=1.-c1
+  DO l=1,ns
+    xx=HUGE(1.)
+    x1=0.0; f1=distrib(x1,c1,c2,sig(l))
+    x2=1.0; f2=distrib(x2,c1,c2,sig(l))
+    DO it=1,maxit
+      x3=0.5*(x1+x2); f3=distrib(x3,c1,c2,sig(l))
+      s=SQRT(f3**2-f1*f2);                 IF(s==0.) EXIT
+      x4=x3+(x3-x1)*(SIGN(1.,f1-f2)*f3/s); IF(ABS(10.*LOG(x4-xx))<=1E-5) EXIT
+      xx=x4; f4=distrib(x4,c1,c2,sig(l));  IF(f4==0.) EXIT
+      IF(SIGN(f3,f4)/=f3) THEN;      x1=x3; f1=f3; x2=xx; f2=f4
+      ELSE IF(SIGN(f1,f4)/=f1) THEN; x2=xx; f2=f4
+      ELSE IF(SIGN(f2,f4)/=f2) THEN; x1=xx; f1=f4
+      ELSE; CALL abort_gcm("ridders",'Algorithm failed (which is odd...')
+      END IF
+      IF(ABS(10.*LOG(ABS(x2-x1)))<=1E-5) EXIT       !--- ERROR ON SIG <= 0.01m            
+    END DO
+    IF(it==maxit+1) WRITE(lunout,'(a,i3)')'WARNING in ridder: failed to converg&
+     &e for level ',l
+    sg(l)=xx
+  END DO
+  sg(1)=1.; sg(ns)=0.
+
+END FUNCTION ridders
+
 END SUBROUTINE disvert
+
+!-------------------------------------------------------------------------------
+
+FUNCTION distrib(x,c1,c2,x0) RESULT(res)
+!
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN) :: x, c1, c2, x0
+  REAL             :: res
+!-------------------------------------------------------------------------------
+  res=c1*x+c2*EXP(1-1/(x**2))-x0
+
+END FUNCTION distrib
+
+