New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
ldfdyn.F90 in branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.
File size: 14.0 KB
Line 
1MODULE ldfdyn
2   !!======================================================================
3   !!                       ***  MODULE  ldfdyn  ***
4   !! Ocean physics:  lateral viscosity coefficient
5   !!=====================================================================
6   !! History :  OPA  ! 1997-07  (G. Madec)  multi dimensional coefficients
7   !!   NEMO     1.0  ! 2002-09  (G. Madec)  F90: Free form and module
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   ldf_dyn_init : initialization, namelist read, and parameters control
12   !!   ldf_dyn_c3d   : 3D eddy viscosity coefficient initialization
13   !!   ldf_dyn_c2d   : 2D eddy viscosity coefficient initialization
14   !!   ldf_dyn_c1d   : 1D eddy viscosity coefficient initialization
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers   
17   USE dom_oce         ! ocean space and time domain
18   USE ldfdyn_oce      ! ocean dynamics lateral physics
19   USE phycst          ! physical constants
20   USE ldfslp          ! ???
21   USE ioipsl
22   USE in_out_manager  ! I/O manager
23   USE lib_mpp         ! distribued memory computing library
24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
25   USE wrk_nemo        ! Memory Allocation
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   ldf_dyn_init   ! called by opa.F90
31
32  INTERFACE ldf_zpf
33     MODULE PROCEDURE ldf_zpf_1d, ldf_zpf_1d_3d, ldf_zpf_3d
34  END INTERFACE
35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38   !!----------------------------------------------------------------------
39   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
40   !! $Id$
41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE ldf_dyn_init
46      !!----------------------------------------------------------------------
47      !!                  ***  ROUTINE ldf_dyn_init  ***
48      !!                   
49      !! ** Purpose :   set the horizontal ocean dynamics physics
50      !!
51      !! ** Method  : 
52      !!      -  default option : ahm = constant coef. = rn_ahm_0 (namelist)
53      !!      - 'key_dynldf_c1d': ahm = F(depth)                     see ldf_dyn_c1d.h90
54      !!      - 'key_dynldf_c2d': ahm = F(latitude,longitude)        see ldf_dyn_c2d.h90
55      !!      - 'key_dynldf_c3d': ahm = F(latitude,longitude,depth)  see ldf_dyn_c3d.h90
56      !!
57      !!      N.B. User defined include files.  By default, 3d and 2d coef.
58      !!      are set to a constant value given in the namelist and the 1d
59      !!      coefficients are initialized to a hyperbolic tangent vertical
60      !!      profile.
61      !!
62      !! Reference :   Madec, G. and M. Imbard, 1996: Climate Dynamics, 12, 381-388.
63      !!----------------------------------------------------------------------
64      INTEGER ::   ioptio         ! ???
65      INTEGER ::   ios            ! Local : output status for namelist read
66      LOGICAL ::   ll_print = .FALSE.    ! Logical flag for printing viscosity coef.
67      !!
68      NAMELIST/namdyn_ldf/ ln_dynldf_lap  , ln_dynldf_bilap,                  &
69         &                 ln_dynldf_level, ln_dynldf_hor  , ln_dynldf_iso,   &
70         &                 rn_ahm_0_lap   , rn_ahmb_0      , rn_ahm_0_blp ,   &
71         &                 rn_cmsmag_1    , rn_cmsmag_2    , rn_cmsh,         &
72         &                 rn_ahm_m_lap   , rn_ahm_m_blp
73
74   !!----------------------------------------------------------------------
75
76      REWIND( numnam_ref )              ! Namelist namdyn_ldf in reference namelist : Lateral physics
77      READ  ( numnam_ref, namdyn_ldf, IOSTAT = ios, ERR = 901)
78901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_ldf in reference namelist', lwp )
79
80      REWIND( numnam_cfg )              ! Namelist namdyn_ldf in configuration namelist : Lateral physics
81      READ  ( numnam_cfg, namdyn_ldf, IOSTAT = ios, ERR = 902 )
82902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_ldf in configuration namelist', lwp )
83      IF(lwm .AND. nprint > 2) WRITE ( numond, namdyn_ldf )
84
85      IF(lwp) THEN                      ! Parameter print
86         WRITE(numout,*)
87         WRITE(numout,*) 'ldf_dyn : lateral momentum physics'
88         WRITE(numout,*) '~~~~~~~'
89         WRITE(numout,*) '   Namelist namdyn_ldf : set lateral mixing parameters'
90         WRITE(numout,*) '      laplacian operator                      ln_dynldf_lap   = ', ln_dynldf_lap
91         WRITE(numout,*) '      bilaplacian operator                    ln_dynldf_bilap = ', ln_dynldf_bilap
92         WRITE(numout,*) '      iso-level                               ln_dynldf_level = ', ln_dynldf_level
93         WRITE(numout,*) '      horizontal (geopotential)               ln_dynldf_hor   = ', ln_dynldf_hor
94         WRITE(numout,*) '      iso-neutral                             ln_dynldf_iso   = ', ln_dynldf_iso
95         WRITE(numout,*) '      horizontal laplacian eddy viscosity     rn_ahm_0_lap    = ', rn_ahm_0_lap
96         WRITE(numout,*) '      background viscosity                    rn_ahmb_0       = ', rn_ahmb_0
97         WRITE(numout,*) '      horizontal bilaplacian eddy viscosity   rn_ahm_0_blp    = ', rn_ahm_0_blp
98         WRITE(numout,*) '      upper limit for laplacian eddy visc     rn_ahm_m_lap    = ', rn_ahm_m_lap
99         WRITE(numout,*) '      upper limit for bilap eddy viscosity    rn_ahm_m_blp    = ', rn_ahm_m_blp
100
101      ENDIF
102
103      ahm0     = rn_ahm_0_lap              ! OLD namelist variables defined from DOCTOR namelist variables
104      ahmb0    = rn_ahmb_0
105      ahm0_blp = rn_ahm_0_blp
106
107      ! ... check of lateral diffusive operator on tracers
108      !           ==> will be done in trazdf module
109
110      ! ... Space variation of eddy coefficients
111      ioptio = 0
112#if defined key_dynldf_c3d
113      IF(lwp) WRITE(numout,*) '   momentum mixing coef. = F( latitude, longitude, depth)'
114      ioptio = ioptio+1
115#endif
116#if defined key_dynldf_c2d
117      IF(lwp) WRITE(numout,*) '   momentum mixing coef. = F( latitude, longitude)'
118      ioptio = ioptio+1
119#endif
120#if defined key_dynldf_c1d
121      IF(lwp) WRITE(numout,*) '   momentum mixing coef. = F( depth )'
122      ioptio = ioptio+1
123      IF( ln_sco ) CALL ctl_stop( 'key_dynldf_c1d cannot be used in s-coordinate (ln_sco)' )
124#endif
125      IF( ioptio == 0 ) THEN
126          IF(lwp) WRITE(numout,*) '   momentum mixing coef. = constant  (default option)'
127        ELSEIF( ioptio > 1 ) THEN
128           CALL ctl_stop( 'use only one of the following keys: key_dynldf_c3d, key_dynldf_c2d, key_dynldf_c1d' )
129      ENDIF
130
131
132      IF( ln_dynldf_bilap ) THEN
133         IF(lwp) WRITE(numout,*) '   biharmonic momentum diffusion'
134         IF( .NOT. ln_dynldf_lap ) ahm0 = ahm0_blp   ! Allow spatially varying coefs, which use ahm0 as input
135         IF( ahm0_blp > 0 .AND. .NOT. lk_esopa )   CALL ctl_stop( 'The horizontal viscosity coef. ahm0 must be negative' )
136      ELSE
137         IF(lwp) WRITE(numout,*) '   harmonic momentum diff. (default)'
138         IF( ahm0 < 0 .AND. .NOT. lk_esopa )   CALL ctl_stop( 'The horizontal viscosity coef. ahm0 must be positive' )
139      ENDIF
140
141
142      ! Lateral eddy viscosity
143      ! ======================
144#if defined key_dynldf_c3d
145      CALL ldf_dyn_c3d( ll_print )   ! ahm = 3D coef. = F( longitude, latitude, depth )
146#elif defined key_dynldf_c2d
147      CALL ldf_dyn_c2d( ll_print )   ! ahm = 1D coef. = F( longitude, latitude )
148#elif defined key_dynldf_c1d
149      CALL ldf_dyn_c1d( ll_print )   ! ahm = 1D coef. = F( depth )
150#else
151                                     ! Constant coefficients
152      IF(lwp) WRITE(numout,*)
153      IF(lwp) WRITE(numout,*) 'inildf: constant eddy viscosity coef. '
154      IF(lwp) WRITE(numout,*) '~~~~~~'
155      IF(lwp) WRITE(numout,*) '        ahm1 = ahm2 = ahm0 =  ',ahm0
156#endif
157     nkahm_smag = 0
158#if defined key_dynldf_smag
159     nkahm_smag = 1
160#endif
161      !
162      IF(lwp .AND. lflush) CALL flush(numout)
163      !
164   END SUBROUTINE ldf_dyn_init
165
166#if defined key_dynldf_c3d
167#   include "ldfdyn_c3d.h90"
168#elif defined key_dynldf_c2d
169#   include "ldfdyn_c2d.h90"
170#elif defined key_dynldf_c1d
171#   include "ldfdyn_c1d.h90"
172#endif
173
174
175   SUBROUTINE ldf_zpf_1d( ld_print, pdam, pwam, pbot, pdep, pah )
176      !!----------------------------------------------------------------------
177      !!                  ***  ROUTINE ldf_zpf  ***
178      !!                   
179      !! ** Purpose :   vertical adimensional profile for eddy coefficient
180      !!
181      !! ** Method  :   1D eddy viscosity coefficients ( depth )
182      !!----------------------------------------------------------------------
183      LOGICAL , INTENT(in   )                 ::   ld_print   ! If true, output arrays on numout
184      REAL(wp), INTENT(in   )                 ::   pdam       ! depth of the inflection point
185      REAL(wp), INTENT(in   )                 ::   pwam       ! width of inflection
186      REAL(wp), INTENT(in   )                 ::   pbot       ! bottom value (0<pbot<= 1)
187      REAL(wp), INTENT(in   ), DIMENSION(jpk) ::   pdep       ! depth of the gridpoint (T, U, V, F)
188      REAL(wp), INTENT(inout), DIMENSION(jpk) ::   pah        ! adimensional vertical profile
189      !!
190      INTEGER  ::   jk           ! dummy loop indices
191      REAL(wp) ::   zm00, zm01, zmhb, zmhs       ! temporary scalars
192      !!----------------------------------------------------------------------
193
194      zm00 = TANH( ( pdam - gdept_1d(1    ) ) / pwam )
195      zm01 = TANH( ( pdam - gdept_1d(jpkm1) ) / pwam )
196      zmhs = zm00 / zm01
197      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
198
199      DO jk = 1, jpk
200         pah(jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam )  )
201      END DO
202
203      IF(lwp .AND. ld_print ) THEN      ! Control print
204         WRITE(numout,*)
205         WRITE(numout,*) '         ahm profile : '
206         WRITE(numout,*)
207         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
208         DO jk = 1, jpk
209            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(jk), pdep(jk)
210         END DO
211      ENDIF
212      !
213      IF(lwp .AND. lflush) CALL flush(numout)
214      !
215   END SUBROUTINE ldf_zpf_1d
216
217
218   SUBROUTINE ldf_zpf_1d_3d( ld_print, pdam, pwam, pbot, pdep, pah )
219      !!----------------------------------------------------------------------
220      !!                  ***  ROUTINE ldf_zpf  ***
221      !!
222      !! ** Purpose :   vertical adimensional profile for eddy coefficient
223      !!
224      !! ** Method  :   1D eddy viscosity coefficients ( depth )
225      !!----------------------------------------------------------------------
226      LOGICAL , INTENT(in   )                         ::   ld_print   ! If true, output arrays on numout
227      REAL(wp), INTENT(in   )                         ::   pdam       ! depth of the inflection point
228      REAL(wp), INTENT(in   )                         ::   pwam       ! width of inflection
229      REAL(wp), INTENT(in   )                         ::   pbot       ! bottom value (0<pbot<= 1)
230      REAL(wp), INTENT(in   ), DIMENSION          (:) ::   pdep       ! depth of the gridpoint (T, U, V, F)
231      REAL(wp), INTENT(inout), DIMENSION      (:,:,:) ::   pah        ! adimensional vertical profile
232      !!
233      INTEGER  ::   jk           ! dummy loop indices
234      REAL(wp) ::   zm00, zm01, zmhb, zmhs, zcf  ! temporary scalars
235      !!----------------------------------------------------------------------
236
237      zm00 = TANH( ( pdam - gdept_1d(1    ) ) / pwam )
238      zm01 = TANH( ( pdam - gdept_1d(jpkm1) ) / pwam )
239      zmhs = zm00 / zm01
240      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
241
242      DO jk = 1, jpk
243         zcf = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam )  )
244         pah(:,:,jk) = zcf
245      END DO
246
247      IF(lwp .AND. ld_print ) THEN      ! Control print
248         WRITE(numout,*)
249         WRITE(numout,*) '         ahm profile : '
250         WRITE(numout,*)
251         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
252         DO jk = 1, jpk
253            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(jk)
254         END DO
255      ENDIF
256      !
257      IF(lwp .AND. lflush) CALL flush(numout)
258      !
259   END SUBROUTINE ldf_zpf_1d_3d
260
261
262   SUBROUTINE ldf_zpf_3d( ld_print, pdam, pwam, pbot, pdep, pah )
263      !!----------------------------------------------------------------------
264      !!                  ***  ROUTINE ldf_zpf  ***
265      !!
266      !! ** Purpose :   vertical adimensional profile for eddy coefficient
267      !!
268      !! ** Method  :   3D for partial step or s-coordinate
269      !!----------------------------------------------------------------------
270      LOGICAL , INTENT(in   )                         ::   ld_print   ! If true, output arrays on numout
271      REAL(wp), INTENT(in   )                         ::   pdam       ! depth of the inflection point
272      REAL(wp), INTENT(in   )                         ::   pwam       ! width of inflection
273      REAL(wp), INTENT(in   )                         ::   pbot       ! bottom value (0<pbot<= 1)
274      REAL(wp), INTENT(in   ), DIMENSION      (:,:,:) ::   pdep       ! dep of the gridpoint (T, U, V, F)
275      REAL(wp), INTENT(inout), DIMENSION      (:,:,:) ::   pah        ! adimensional vertical profile
276      !!
277      INTEGER  ::   jk           ! dummy loop indices
278      REAL(wp) ::   zm00, zm01, zmhb, zmhs       ! temporary scalars
279      !!----------------------------------------------------------------------
280
281      zm00 = TANH( ( pdam - gdept_1d(1    ) ) / pwam )   
282      zm01 = TANH( ( pdam - gdept_1d(jpkm1) ) / pwam )
283      zmhs = zm00 / zm01
284      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
285
286      DO jk = 1, jpk
287         pah(:,:,jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(:,:,jk) ) / pwam )  )
288      END DO
289
290      IF(lwp .AND. ld_print ) THEN      ! Control print
291         WRITE(numout,*)
292         WRITE(numout,*) '         ahm profile : '
293         WRITE(numout,*)
294         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
295         DO jk = 1, jpk
296            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(1,1,jk)
297         END DO
298      ENDIF
299      !
300      IF(lwp .AND. lflush) CALL flush(numout)
301      !
302   END SUBROUTINE ldf_zpf_3d
303
304   !!======================================================================
305END MODULE ldfdyn
Note: See TracBrowser for help on using the repository browser.