Changeset 2379 for branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

Timestamp:

2010-11-12T15:51:18+01:00 (13 years ago)

Author:

acc

Message:

nemo_v3_3beta. First stage of adding ORCA1 support. See ticket #758. Added optional double tanh function for vertical grid

File:

• branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 (modified) (8 diffs)

branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

@@ -160 +160 @@
       REAL(wp) ::   zacr, zdzmin, zhmax    ! par_CONFIG_Rxx.h90
       REAL(wp) ::   zrefdep                ! depth of the reference level (~10m)
+      REAL(wp) ::   za2, zkth2, zacr2      ! Values for optional double tanh function set from parameters 
       !!----------------------------------------------------------------------
 
@@ -166 +167 @@
        zkth = ppkth       ;   zacr = ppacr
        zdzmin = ppdzmin   ;   zhmax = pphmax
+       zkth2 = ppkth2     ;   zacr2 = ppacr2   ! optional (ldbletanh=T) double tanh parameters
 
       ! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
@@ -181 +183 @@
       ELSE
          za1 = ppa1 ;       za0 = ppa0 ;          zsur = ppsur
+         za2 = ppa2                            ! optional (ldbletanh=T) double tanh parameter
       ENDIF
 
@@ -203 +206 @@
             WRITE(numout,*) '                 zkth = ', zkth
             WRITE(numout,*) '                 zacr = ', zacr
+            IF( ldbletanh ) THEN
+               WRITE(numout,*) ' (Double tanh    za2  = ', za2
+               WRITE(numout,*) '  parameters)    zkth2= ', zkth2
+               WRITE(numout,*) '                 zacr2= ', zacr2
+            ENDIF
          ENDIF
       ENDIF
@@ -220 +228 @@
          END DO
       ELSE                                ! Madec & Imbard 1996 function
-         DO jk = 1, jpk
-            zw = FLOAT( jk )
-            zt = FLOAT( jk ) + 0.5
-            gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  )
-            gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  )
-            e3w_0  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )
-            e3t_0  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )
-         END DO
+         IF( .NOT. ldbletanh ) THEN
+            DO jk = 1, jpk
+               zw = FLOAT( jk )
+               zt = FLOAT( jk ) + 0.5
+               gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  )
+               gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  )
+               e3w_0  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )
+               e3t_0  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )
+            END DO
+         ELSE
+            DO jk = 1, jpk
+               zw = FLOAT( jk )
+               zt = FLOAT( jk ) + 0.5
+               ! Double tanh function
+               gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )    &
+     &                                         + za2 * zacr2* LOG ( COSH( (zw-zkth2)/zacr2 ) )  )
+               gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )    &
+     &                                         + za2 * zacr2* LOG ( COSH( (zt-zkth2)/zacr2 ) )  )
+               e3w_0  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )    &
+     &                                         + za2        * TANH(       (zw-zkth2)/zacr2   )
+               e3t_0  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )    &
+     &                                         + za2        * TANH(       (zt-zkth2)/zacr2   )
+            END DO
+         ENDIF
          gdepw_0(1) = 0.e0                     ! force first w-level to be exactly at zero
       ENDIF
 
 !!gm BUG in s-coordinate this does not work!
-      ! deepest/shallowest W level Above/Bellow ~10m
+      ! deepest/shallowest W level Above/Below ~10m
       zrefdep = 10. - ( 0.1*MINVAL(e3w_0) )                          ! ref. depth with tolerance (10% of minimum layer thickness)
-      nlb10 = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 )   ! shallowest W level Bellow ~10m
-      nla10 = nlb10 - 1                                              ! deepest    W level Above  ~10m
+      nlb10 = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 )   ! shallowest W level Below ~10m
+      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
 !!gm end bug
 
@@ -746 +770 @@
       INTEGER  ::   ji, jj, jk       ! dummy loop indices
       INTEGER  ::   ik, it           ! temporary integers
+      INTEGER  ::   nlbelow          ! temporary integer
       LOGICAL  ::   ll_print         ! Allow  control print for debugging
       REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
@@ -751 +776 @@
       REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth
       REAL(wp) ::   zdiff            ! temporary scalar
+      REAL(wp) ::   zrefdep          ! temporary scalar
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zprt   ! 3D workspace
       !!---------------------------------------------------------------------
@@ -772 +798 @@
       ! ===================
       zmax = gdepw_0(jpk) + e3t_0(jpk)     ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) )
-      zmin = gdepw_0(4)                    ! minimum depth = 3 levels
+
+      zrefdep = 25.
+      nlbelow = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 )   ! shallowest W level Below zrefdep
+      IF(lwp) write(numout,*) 'Minimum depth level selected: ', nlbelow
+      zmin = gdepw_0(nlbelow)              ! minimum depth = nlbelow-1 levels
 
       mbathy(:,:) = jpkm1                  ! initialize mbathy to the maximum ocean level available