Changeset 606

Timestamp:

10/25/17 17:00:32 (7 years ago)

Author:

dubos

Message:

trunk : new disvert strato_custom

Location:

codes/icosagcm/trunk/src/vertical

Files:

• disvert.f90 (modified) (2 diffs)
• disvert_strato.f90 (modified) (3 diffs)

codes/icosagcm/trunk/src/vertical/disvert.f90

@@ -23 +23 @@
   SUBROUTINE init_disvert
   USE disvert_std_mod, ONLY: ap_std=>ap, bp_std=>bp, presnivs_std=>presnivs, init_disvert_std=>init_disvert
-  USE disvert_strato_mod, ONLY: ap_strato=>ap, bp_strato=>bp, presnivs_strato=>presnivs, init_disvert_strato=>init_disvert
+  USE disvert_strato_mod, ONLY: ap_strato=>ap, bp_strato=>bp, presnivs_strato=>presnivs, init_disvert_strato=>init_disvert, init_disvert_strato_custom
   USE disvert_apbp_mod, ONLY: ap_apbp=>ap, bp_apbp=>bp, presnivs_apbp=>presnivs, init_disvert_apbp=>init_disvert
   USE disvert_ncar_mod, ONLY: ap_ncar=>ap, bp_ncar=>bp, presnivs_ncar=>presnivs, init_disvert_ncar=>init_disvert
@@ -69 +69 @@
         
         CALL init_disvert_strato
+        ap=>ap_strato
+        bp=>bp_strato
+        presnivs=>presnivs_strato
+
+      CASE('strato_custom')
+        
+        CALL init_disvert_strato_custom
         ap=>ap_strato
         bp=>bp_strato

codes/icosagcm/trunk/src/vertical/disvert_strato.f90

@@ -1 +1 @@
 MODULE disvert_strato_mod
   USE icosa
+  USE mpipara
+  USE getin_mod
+  USE earth_const
+  IMPLICIT NONE
+  PRIVATE
+
   REAL(rstd), SAVE, ALLOCATABLE,TARGET :: ap(:)
 !$OMP THREADPRIVATE(ap)
@@ -8 +14 @@
 !$OMP THREADPRIVATE(presnivs)
 
+  PUBLIC :: ap,bp,presnivs, init_disvert, init_disvert_strato_custom
+
 CONTAINS
-
+  
   SUBROUTINE init_disvert
-  USE icosa
-  USE getin_mod
-  IMPLICIT NONE
-  REAL(rstd) :: dsigmin
+    REAL(rstd) :: dsigmin, snorm, x, dsig(llm),  sig(llm+1)
+    INTEGER :: l
+
     ALLOCATE(ap(llm+1))
     ALLOCATE(bp(llm+1))
     ALLOCATE(presnivs(llm))
-
+    
     dsigmin=1.0  ! Should be 0.3 for CMIP5
     CALL getin('disvert_dsigmin', dsigmin) 
-    CALL disvert(dsigmin,ap,bp,presnivs)
-
-  END SUBROUTINE init_disvert  
-
-
-  SUBROUTINE disvert(dsigmin,ap,bp,presnivs)
-  USE icosa
-  USE mpipara
-  USE earth_const
-  IMPLICIT NONE
-  REAL(rstd) :: dsigmin
-  REAL(rstd),INTENT(OUT) :: ap(:)
-  REAL(rstd),INTENT(OUT) :: bp(:)
-  REAL(rstd),INTENT(OUT) :: presnivs(:)
-  
-  REAL(rstd) :: dsig(llm)
-  REAL(rstd) :: sig(llm+1)
-  REAL(rstd) :: snorm
-  REAL(rstd) :: x
-  INTEGER :: l
-  
+    
     snorm  = 0.
     DO l = 1, llm
@@ -84 +71 @@
     ENDIF
   
-  END SUBROUTINE disvert
+  END SUBROUTINE init_disvert
+  
+  SUBROUTINE init_disvert_strato_custom
+    REAL(rstd) :: vert_scale,vert_dzmin,vert_dzlow,vert_z0low,vert_dzmid,vert_z0mid,vert_h_mid,vert_dzhig,vert_z0hig,vert_h_hig
+    REAL(rstd) :: z, snorm, dsig(llm), sig(llm+1), sig0(llm+1), zz(llm+1)
+    INTEGER :: l
+    ALLOCATE(ap(llm+1))
+    ALLOCATE(bp(llm+1))
+    ALLOCATE(presnivs(llm))
+!===================================================================
+! David Cugnet, François Lott, Lionel Guez, Ehouoarn Millour, Fredho
+! 2014/05
+! custumize strato distribution
+! Al the parameter are given in km assuming a given scalehigh
+    vert_scale=7.     ! scale hight
+    vert_dzmin=0.02   ! width of first layer
+    vert_dzlow=1.     ! dz in the low atmosphere
+    vert_z0low=8.     ! height at which resolution recches dzlow
+    vert_dzmid=3.     ! dz in the mid atmsophere
+    vert_z0mid=70.    ! height at which resolution recches dzmid
+    vert_h_mid=20.    ! width of the transition
+    vert_dzhig=11.    ! dz in the high atmsophere
+    vert_z0hig=80.    ! height at which resolution recches dz
+    vert_h_hig=20.    ! width of the transition
+!===================================================================
+
+    call getin('vert_scale',vert_scale)
+    call getin('vert_dzmin',vert_dzmin)
+    call getin('vert_dzlow',vert_dzlow)
+    call getin('vert_z0low',vert_z0low)
+    CALL getin('vert_dzmid',vert_dzmid)
+    CALL getin('vert_z0mid',vert_z0mid)
+    call getin('vert_h_mid',vert_h_mid)
+    call getin('vert_dzhig',vert_dzhig)
+    call getin('vert_z0hig',vert_z0hig)
+    call getin('vert_h_hig',vert_h_hig)
+  
+!!!!!    scaleheight=vert_scale ! for consistency with further computations
+    ! Build geometrical distribution
+    zz(1)=0.
+    DO l=1,llm
+       z=zz(l)
+       zz(l+1)=zz(l)+vert_dzmin+vert_dzlow*TANH(z/vert_z0low)+                &
+                  (vert_dzmid-vert_dzlow)* &
+                       (TANH((z-vert_z0mid)/vert_h_mid)-TANH((-vert_z0mid)/vert_h_mid)) &
+                 +(vert_dzhig-vert_dzmid-vert_dzlow)*                                  &
+                       (TANH((z-vert_z0hig)/vert_h_hig)-TANH((-vert_z0hig)/vert_h_hig))
+    ENDDO
+
+
+!===================================================================
+! Comment added Fredho 2014/05/18, Saint-Louis, Senegal
+! From approximate z to ap, bp, so that p=ap+bp*p0 and p/p0=exp(-z/H)
+! sig0 is p/p0
+! sig is an intermediate distribution introduce to estimate bp
+! 1.   sig0=exp(-z/H)
+! 2.   inversion of sig0=(1-pa/p0)*sig+(1-pa/p0)*exp(1-1/sig**2)
+! 3.   bp=exp(1-1/sig**2)
+! 4.   ap deduced from  the combination of 2 and 3 so that sig0=ap/p0+bp
+!===================================================================
+
+    sig0(llm+1)=0.
+    DO l = 1, llm
+       sig0(l)=EXP(-zz(l)/vert_scale)
+!!!     sig(l)=racinesig(sig0(l))
+    ENDDO
+    sig=racinesig(sig0)
+    
+    bp(llm+1) =   0.
+    DO l = 2, llm
+       bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
+       ap(l) = pa * ( sig(l) - bp(l) )
+    ENDDO
+    bp(1)=1.
+    !!    ap=pa*(sig-bp)
+    ap(1)=0.
+    ap(llm+1) = pa * ( sig(llm+1) - bp(llm+1) )
+    DO l = 1, llm
+       presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
+    ENDDO
+    
+    ! tell the world about it
+    IF (is_mpi_root) THEN
+       !$OMP MASTER
+       WRITE(*,*) "ap()=",ap
+       WRITE(*,*) "bp()=",bp
+       WRITE(*,*) "Approximative mid-layer pressure, assuming a surface pressure preff=",preff," Pa"
+       WRITE(*,*) "and approximative mid-layer height, assuming an atmospheric scale height of ",scale_height/1000," (km)"
+       DO l=1,llm
+          WRITE(*,*) 'PRESNIVS(',l,')=',presnivs(l),'  Z ~ ',log(preff/presnivs(l))*scale_height/1000,       &
+               ' DZ ~ ',scale_height/1000*log((ap(l)+bp(l)*preff)/ max(ap(l+1)+bp(l+1)*preff,1.e-10))
+       ENDDO
+       !$OMP END MASTER
+    ENDIF
+    
+  CONTAINS
+    
+    FUNCTION racinesig(sig) RESULT(sg)
+      !
+      !-------------------------------------------------------------------------------
+      ! Fredho 2014/05/18
+      ! Purpose: Search for s solving (Pa/Preff)*sg+(1-Pa/Preff)*EXP(1-1./sg**2)=s
+      ! Notes:   Uses Newton Raphson search
+      !-------------------------------------------------------------------------------
+      ! Arguments:
+      REAL, INTENT(IN)  :: sig(:)
+      REAL(rstd)        :: sg(SIZE(sig))
+      !-------------------------------------------------------------------------------
+      ! Local variables:
+      INTEGER      :: it, ns, maxit
+      REAL(rstd)   :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
+      !-------------------------------------------------------------------------------
+      ns=SIZE(sig); maxit=100
+      c1=Pa/Preff; c2=1.-c1
+      DO l=2,ns-1
+         sg(l)=sig(l)
+         DO it=1,maxit
+            f1=exp(1-1./sg(l)**2)*(1.-c1)
+            sg(l)=sg(l)-(c1*sg(l)+f1-sig(l))/(c1+2*f1*sg(l)**(-3))
+         ENDDO
+         !   print*,'SSSSIG ',sig(l),sg(l),c1*sg(l)+exp(1-1./sg(l)**2)*(1.-c1)
+      ENDDO
+      sg(1)=1.; sg(ns)=0.
+      
+    END FUNCTION racinesig
+    
+  END SUBROUTINE init_disvert_strato_custom
   
 END MODULE disvert_strato_mod