Changeset 1839 for branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF

Timestamp:

2010-04-14T15:28:04+02:00 (14 years ago)

Author:

edblockley

Message:

1st commit for the mixed_ldfdyn branch
see ticket:659

Location:

branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF

Files:

• ldfdyn.F90 (modified) (3 diffs)
• ldfdyn_c1d.h90 (modified) (5 diffs)
• ldfdyn_c2d.h90 (modified) (2 diffs)
• ldfdyn_c3d.h90 (modified) (1 diff)
• ldfdyn_oce.F90 (modified) (1 diff)
• ldfdyn_substitute.h90 (modified) (2 diffs)

branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn.F90

@@ -67 +67 @@
       NAMELIST/namdyn_ldf/ ln_dynldf_lap  , ln_dynldf_bilap,                  &
          &                 ln_dynldf_level, ln_dynldf_hor  , ln_dynldf_iso,   &
-         &                 rn_ahm_0       , rn_ahmb_0
+         &                 rn_ahm_0       , rn_ahmb_0      , rn_ahm_0_blp
       !!----------------------------------------------------------------------
 
@@ -78 +78 @@
          WRITE(numout,*) '~~~~~~~'
          WRITE(numout,*) '   Namelist nam_dynldf : set lateral mixing parameters'
-         WRITE(numout,*) '      laplacian operator          ln_dynldf_lap   = ', ln_dynldf_lap
-         WRITE(numout,*) '      bilaplacian operator        ln_dynldf_bilap = ', ln_dynldf_bilap
-         WRITE(numout,*) '      iso-level                   ln_dynldf_level = ', ln_dynldf_level
-         WRITE(numout,*) '      horizontal (geopotential)   ln_dynldf_hor   = ', ln_dynldf_hor
-         WRITE(numout,*) '      iso-neutral                 ln_dynldf_iso   = ', ln_dynldf_iso
-         WRITE(numout,*) '      horizontal eddy viscosity   rn_ahm_0        = ', rn_ahm_0
-         WRITE(numout,*) '      background viscosity        rn_ahmb_0       = ', rn_ahmb_0
-      ENDIF
-
-      ahm0  = rn_ahm_0                  ! OLD namelist variables defined from DOCTOR namelist variables
-      ahmb0 = rn_ahmb_0
+         WRITE(numout,*) '      laplacian operator                      ln_dynldf_lap   = ', ln_dynldf_lap
+         WRITE(numout,*) '      bilaplacian operator                    ln_dynldf_bilap = ', ln_dynldf_bilap
+         WRITE(numout,*) '      iso-level                               ln_dynldf_level = ', ln_dynldf_level
+         WRITE(numout,*) '      horizontal (geopotential)               ln_dynldf_hor   = ', ln_dynldf_hor
+         WRITE(numout,*) '      iso-neutral                             ln_dynldf_iso   = ', ln_dynldf_iso
+         WRITE(numout,*) '      horizontal laplacian eddy viscosity     rn_ahm_0        = ', rn_ahm_0
+         WRITE(numout,*) '      background viscosity                    rn_ahmb_0       = ', rn_ahmb_0
+         WRITE(numout,*) '      horizontal bilaplacian eddy viscosity   rn_ahm_0        = ', rn_ahm_0
+
+      ENDIF
+
+      ahm0     = rn_ahm_0                  ! OLD namelist variables defined from DOCTOR namelist variables
+      ahmb0    = rn_ahmb_0
+      ahm0_blp = rn_ahm_0_blp
 
       ! ... check of lateral diffusive operator on tracers
@@ -117 +120 @@
       IF( ln_dynldf_bilap ) THEN
          IF(lwp) WRITE(numout,*) '   biharmonic momentum diffusion'
-         IF( ahm0 > 0 .AND. .NOT. lk_esopa )   CALL ctl_stop( 'The horizontal viscosity coef. ahm0 must be negative' )
+         IF( ahm0_blp == 0.0 ) ahm0_blp = ahm0       ! Old namelist method: bilap specified with ahm0
+         IF( .NOT. ln_dynldf_lap ) ahm0 = ahm0_blp   ! Allow spatially varying coefs, which use ahm0 as input
+         IF( ahm0_blp > 0 .AND. .NOT. lk_esopa )   CALL ctl_stop( 'The horizontal viscosity coef. ahm0 must be negative' )
       ELSE
          IF(lwp) WRITE(numout,*) '   harmonic momentum diff. (default)'

branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_c1d.h90

@@ -17 +17 @@
       !! ** Method  :   1D eddy viscosity coefficients ( depth )
       !!      ahm3, ahm4 never used
-      !!      biharmonic or harmonic operator : ahm1=ahm2 defined at T-level
+      !!      harmonic   operator : ahm1 defined at T-level
+      !!      biharmonic operator : ahm2 defined at T-level
       !!      isopycnal or geopotential harmonic operator
       !!                                      : ahm1 defined at T-level
@@ -27 +28 @@
 
       !! * Local variables
-      REAL(wp) ::   zdam, zwam, zm00, zm01, zmhf, zmhs
+      REAL(wp) ::    zdam,  zwam,  zm00,  zm01,  zmhf,  zmhs
+      REAL(wp) ::   zdam2, zwam2, zm200, zm201, zmh2f, zmh2s
       REAL(wp) ::   zahmf, zahms
       !!----------------------------------------------------------------------
@@ -36 +38 @@
       IF(lwp) WRITE(numout,*)
 
-      ! Set ahm1=ahm2   (always at t-level)
+      ! Set ahm1 for laplacian     (always at t-level)
       ! =============
       ! (USER: modify ahm1 following your desiderata)
@@ -53 +55 @@
       zmhs = zahms-zmhf * zm00
 
-      ! set ahm1=ahm2 at T-level
+      ! Set ahm2 for bilaplacian   (always at t-level)
+      ! =============
+      ! (USER: modify ahm2 following your desiderata)
+
+      ! initialization of the profile
+      ! ahms, ahmf: surface and bottom values
+      zahm2s = ahm0_blp
+      zahm2f = ahm0_blp/4.
+      ! zdam, zwam: depth of the inflection pt and width of inflection
+      zdam2  = -300.
+      zwam2  =  300.
+      ! computation coefficients
+      zm200 = TANH( (0-zdam2)/zwam2 )
+      zm201 = TANH( (-fsdept(1,1,jpk)-zdam2)/zwam2 )
+      zmh2f = (zahm2s-zahm2f)/(zm200-zm201)
+      zmh2s = zahm2s-zmh2f * zm00
+
+
+      ! set ahm1 and ahm2 at T-level
       DO jk = 1, jpk
-         ahm1(jk) = zmhs + zmhf * TANH( (-fsdept(1,1,jk)-zdam) / zwam )
-         ahm2(jk) = ahm1(jk)
+         ahm1(jk) = zmhs  + zmhf  * TANH( (-fsdept(1,1,jk)-zdam ) / zwam  )
+         ahm2(jk) = zmh2s + zmh2f * TANH( (-fsdept(1,1,jk)-zdam2) / zwam2 )
       END DO
 
@@ -62 +82 @@
       IF(lwp .AND. ld_print ) THEN
          WRITE(numout,*)
-         WRITE(numout,*) '         ahm profile : '
+         WRITE(numout,*) '         ahm profile (laplacian): '
          WRITE(numout,*)
          WRITE(numout,9100)
          DO jk = 1, jpk
             WRITE(numout,9110) jk, ahm1(jk), fsdept(1,1,jk)
+         END DO
+         WRITE(numout,*)
+         WRITE(numout,*) '         ahm profile (bilaplacian): '
+         WRITE(numout,*)
+         WRITE(numout,9100)
+         DO jk = 1, jpk
+            WRITE(numout,9110) jk, ahm2(jk), fsdept(1,1,jk)
          END DO
       ENDIF

branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90

@@ -27 +27 @@
       !!           iso-model level : ahm3, ahm4 not used
       !!
-      !!       biharmonic operator : ahm1 is defined at u-point
-      !!                             ahm2 is defined at v-point
-      !!                           : ahm3, ahm4 not used
+      !!       biharmonic operator : ahm3 is defined at u-point
+      !!                             ahm4 is defined at v-point
+      !!                           : ahm1, ahm2 not used
       !!
       !!----------------------------------------------------------------------
@@ -102 +102 @@
          IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
 
-         za00 = ahm0 / ( zd_max * zd_max * zd_max )
+         za00 = ahm0_blp / ( zd_max * zd_max * zd_max )
          DO jj = 1, jpj
             DO ji = 1, jpi

branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90

@@ -122 +122 @@
          IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
 
-         za00 = ahm0 / ( zd_max * zd_max * zd_max )
+         za00 = ahm0_blp / ( zd_max * zd_max * zd_max )
          DO jj = 1, jpj
             DO ji = 1, jpi

branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_oce.F90

@@ -17 +17 @@
    LOGICAL , PUBLIC ::   ln_dynldf_hor   = .TRUE.      !: horizontal (geopotential) direction
    LOGICAL , PUBLIC ::   ln_dynldf_iso   = .FALSE.     !: iso-neutral direction
-   REAL(wp), PUBLIC ::   rn_ahm_0        = 40000._wp   !: lateral eddy viscosity (m2/s)
-   REAL(wp), PUBLIC ::   rn_ahmb_0       =     0._wp   !: lateral background eddy viscosity (m2/s)
+   REAL(wp), PUBLIC ::   rn_ahm_0        = 40000._wp   !: lateral laplacian eddy viscosity (m2/s)
+   REAL(wp), PUBLIC ::   rn_ahmb_0       =     0._wp   !: lateral laplacian background eddy viscosity (m2/s)
+   REAL(wp), PUBLIC ::   rn_ahm_0_blp    =     0._wp   !: lateral bilaplacian eddy viscosity (m4/s)
 
-   REAL(wp), PUBLIC ::   ahm0, ahmb0                 ! OLD namelist names
+   REAL(wp), PUBLIC ::   ahm0, ahmb0, ahm0_blp       ! OLD namelist names
 
 #if defined key_dynldf_c3d

branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_substitute.h90

@@ -10 +10 @@
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    !!----------------------------------------------------------------------
+   !!
+   !! fsahmt, fsahmf - used for laplaian operator only
+   !! fsahmu, fsahmv - used for bilaplacian operator only
+   !!
 #if defined key_dynldf_c3d
 !   ' key_dynldf_c3d' :                  3D coefficient
@@ -26 +30 @@
 #  define   fsahmt(i,j,k)    ahm1(k)
 #  define   fsahmf(i,j,k)    ahm1(k)
-#  define   fsahmu(i,j,k)    ahm1(k)
-#  define   fsahmv(i,j,k)    ahm1(k)
+#  define   fsahmu(i,j,k)    ahm2(k)
+#  define   fsahmv(i,j,k)    ahm2(k)
 #else
 !   default option :               Constant coefficient
 #  define   fsahmt(i,j,k)   ahm0
 #  define   fsahmf(i,j,k)   ahm0
-#  define   fsahmu(i,j,k)   ahm0
-#  define   fsahmv(i,j,k)   ahm0
+#  define   fsahmu(i,j,k)   ahm0_blp
+#  define   fsahmv(i,j,k)   ahm0_blp
 #endif