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 trunk/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 4624

Last change on this file since 4624 was 4624, checked in by acc, 10 years ago

#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.
  • Property svn:keywords set to Id
File size: 13.8 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) 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      !
163   END SUBROUTINE ldf_dyn_init
164
165#if defined key_dynldf_c3d
166#   include "ldfdyn_c3d.h90"
167#elif defined key_dynldf_c2d
168#   include "ldfdyn_c2d.h90"
169#elif defined key_dynldf_c1d
170#   include "ldfdyn_c1d.h90"
171#endif
172
173
174   SUBROUTINE ldf_zpf_1d( ld_print, pdam, pwam, pbot, pdep, pah )
175      !!----------------------------------------------------------------------
176      !!                  ***  ROUTINE ldf_zpf  ***
177      !!                   
178      !! ** Purpose :   vertical adimensional profile for eddy coefficient
179      !!
180      !! ** Method  :   1D eddy viscosity coefficients ( depth )
181      !!----------------------------------------------------------------------
182      LOGICAL , INTENT(in   )                 ::   ld_print   ! If true, output arrays on numout
183      REAL(wp), INTENT(in   )                 ::   pdam       ! depth of the inflection point
184      REAL(wp), INTENT(in   )                 ::   pwam       ! width of inflection
185      REAL(wp), INTENT(in   )                 ::   pbot       ! bottom value (0<pbot<= 1)
186      REAL(wp), INTENT(in   ), DIMENSION(jpk) ::   pdep       ! depth of the gridpoint (T, U, V, F)
187      REAL(wp), INTENT(inout), DIMENSION(jpk) ::   pah        ! adimensional vertical profile
188      !!
189      INTEGER  ::   jk           ! dummy loop indices
190      REAL(wp) ::   zm00, zm01, zmhb, zmhs       ! temporary scalars
191      !!----------------------------------------------------------------------
192
193      zm00 = TANH( ( pdam - gdept_1d(1    ) ) / pwam )
194      zm01 = TANH( ( pdam - gdept_1d(jpkm1) ) / pwam )
195      zmhs = zm00 / zm01
196      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
197
198      DO jk = 1, jpk
199         pah(jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam )  )
200      END DO
201
202      IF(lwp .AND. ld_print ) THEN      ! Control print
203         WRITE(numout,*)
204         WRITE(numout,*) '         ahm profile : '
205         WRITE(numout,*)
206         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
207         DO jk = 1, jpk
208            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(jk), pdep(jk)
209         END DO
210      ENDIF
211      !
212   END SUBROUTINE ldf_zpf_1d
213
214
215   SUBROUTINE ldf_zpf_1d_3d( ld_print, pdam, pwam, pbot, pdep, pah )
216      !!----------------------------------------------------------------------
217      !!                  ***  ROUTINE ldf_zpf  ***
218      !!
219      !! ** Purpose :   vertical adimensional profile for eddy coefficient
220      !!
221      !! ** Method  :   1D eddy viscosity coefficients ( depth )
222      !!----------------------------------------------------------------------
223      LOGICAL , INTENT(in   )                         ::   ld_print   ! If true, output arrays on numout
224      REAL(wp), INTENT(in   )                         ::   pdam       ! depth of the inflection point
225      REAL(wp), INTENT(in   )                         ::   pwam       ! width of inflection
226      REAL(wp), INTENT(in   )                         ::   pbot       ! bottom value (0<pbot<= 1)
227      REAL(wp), INTENT(in   ), DIMENSION          (:) ::   pdep       ! depth of the gridpoint (T, U, V, F)
228      REAL(wp), INTENT(inout), DIMENSION      (:,:,:) ::   pah        ! adimensional vertical profile
229      !!
230      INTEGER  ::   jk           ! dummy loop indices
231      REAL(wp) ::   zm00, zm01, zmhb, zmhs, zcf  ! temporary scalars
232      !!----------------------------------------------------------------------
233
234      zm00 = TANH( ( pdam - gdept_1d(1    ) ) / pwam )
235      zm01 = TANH( ( pdam - gdept_1d(jpkm1) ) / pwam )
236      zmhs = zm00 / zm01
237      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
238
239      DO jk = 1, jpk
240         zcf = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam )  )
241         pah(:,:,jk) = zcf
242      END DO
243
244      IF(lwp .AND. ld_print ) THEN      ! Control print
245         WRITE(numout,*)
246         WRITE(numout,*) '         ahm profile : '
247         WRITE(numout,*)
248         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
249         DO jk = 1, jpk
250            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(jk)
251         END DO
252      ENDIF
253      !
254   END SUBROUTINE ldf_zpf_1d_3d
255
256
257   SUBROUTINE ldf_zpf_3d( ld_print, pdam, pwam, pbot, pdep, pah )
258      !!----------------------------------------------------------------------
259      !!                  ***  ROUTINE ldf_zpf  ***
260      !!
261      !! ** Purpose :   vertical adimensional profile for eddy coefficient
262      !!
263      !! ** Method  :   3D for partial step or s-coordinate
264      !!----------------------------------------------------------------------
265      LOGICAL , INTENT(in   )                         ::   ld_print   ! If true, output arrays on numout
266      REAL(wp), INTENT(in   )                         ::   pdam       ! depth of the inflection point
267      REAL(wp), INTENT(in   )                         ::   pwam       ! width of inflection
268      REAL(wp), INTENT(in   )                         ::   pbot       ! bottom value (0<pbot<= 1)
269      REAL(wp), INTENT(in   ), DIMENSION      (:,:,:) ::   pdep       ! dep of the gridpoint (T, U, V, F)
270      REAL(wp), INTENT(inout), DIMENSION      (:,:,:) ::   pah        ! adimensional vertical profile
271      !!
272      INTEGER  ::   jk           ! dummy loop indices
273      REAL(wp) ::   zm00, zm01, zmhb, zmhs       ! temporary scalars
274      !!----------------------------------------------------------------------
275
276      zm00 = TANH( ( pdam - gdept_1d(1    ) ) / pwam )   
277      zm01 = TANH( ( pdam - gdept_1d(jpkm1) ) / pwam )
278      zmhs = zm00 / zm01
279      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
280
281      DO jk = 1, jpk
282         pah(:,:,jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(:,:,jk) ) / pwam )  )
283      END DO
284
285      IF(lwp .AND. ld_print ) THEN      ! Control print
286         WRITE(numout,*)
287         WRITE(numout,*) '         ahm profile : '
288         WRITE(numout,*)
289         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
290         DO jk = 1, jpk
291            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(1,1,jk)
292         END DO
293      ENDIF
294      !
295   END SUBROUTINE ldf_zpf_3d
296
297   !!======================================================================
298END MODULE ldfdyn
Note: See TracBrowser for help on using the repository browser.