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

source: trunk/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 461

Last change on this file since 461 was 461, checked in by opalod, 18 years ago

nemo_v1_update_052:RB: update lateral diffusion computation following the reorganization of both dynamics and tracers
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 13.2 KB
Line 
1MODULE ldfdyn
2   !!======================================================================
3   !!                       ***  MODULE  ldfdyn  ***
4   !! Ocean physics:  lateral viscosity coefficient
5   !!=====================================================================
6
7   !!----------------------------------------------------------------------
8   !!   ldf_dyn_init : initialization, namelist read, and parameters control
9   !!   ldf_dyn_c3d   : 3D eddy viscosity coefficient initialization
10   !!   ldf_dyn_c2d   : 2D eddy viscosity coefficient initialization
11   !!   ldf_dyn_c1d   : 1D eddy viscosity coefficient initialization
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE oce             ! ocean dynamics and tracers   
15   USE dom_oce         ! ocean space and time domain
16   USE ldfdyn_oce      ! ocean dynamics lateral physics
17   USE phycst          ! physical constants
18   USE ldfslp          ! ???
19   USE in_out_manager  ! I/O manager
20   USE lib_mpp         ! distribued memory computing library
21   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
22
23   IMPLICIT NONE
24   PRIVATE
25
26   !! *  Routine accessibility
27   PUBLIC ldf_dyn_init   ! called by opa.F90
28
29  INTERFACE ldf_zpf
30     MODULE PROCEDURE ldf_zpf_1d, ldf_zpf_1d_3d, ldf_zpf_3d
31  END INTERFACE
32
33   !! * Substitutions
34#  include "domzgr_substitute.h90"
35   !!----------------------------------------------------------------------
36   !!   OPA 9.0 , LOCEAN-IPSL (2005)
37   !! $Header$
38   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
39   !!----------------------------------------------------------------------
40
41CONTAINS
42
43   SUBROUTINE ldf_dyn_init
44      !!----------------------------------------------------------------------
45      !!                  ***  ROUTINE ldf_dyn_init  ***
46      !!                   
47      !! ** Purpose :   set the horizontal ocean dynamics physics
48      !!
49      !! ** Method  : 
50      !!      Eddy viscosity coefficients:
51      !!         default option   : constant coef. ahm0 (namelist)
52      !!        'key_dynldf_c1d': depth dependent coef. defined in
53      !!                        in ldf_dyn_c1d routine
54      !!        'key_dynldf_c2d': latitude and longitude dependent coef.
55      !!                        defined in ldf_dyn_c2d routine
56      !!        'key_dynldf_c3d': latitude, longitude, depth dependent coef.
57      !!                        defined in ldf_dyn_c3d routine
58      !!      N.B. User defined include files.  By default, 3d and 2d coef.
59      !!      are set to a constant value given in the namelist and the 1d
60      !!      coefficients are initialized to a hyperbolic tangent vertical
61      !!      profile.
62      !!
63      !! Reference :
64      !!      Madec, G. and M. Imbard, 1996, A global ocean mesh to overcome
65      !!      the North Pole singularity, Climate Dynamics, 12, 381-388.
66      !!
67      !! History :
68      !!        !  07-97  (G. Madec)  from inimix.F split in 2 routines
69      !!        !  08-97  (G. Madec)  multi dimensional coefficients
70      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
71      !!----------------------------------------------------------------------
72      !! * Modules used
73      USE ioipsl
74
75      !! * Local declarations
76      INTEGER ::   ioptio         ! ???
77      LOGICAL :: ll_print = .FALSE.    ! Logical flag for printing viscosity coef.
78
79       
80      NAMELIST/nam_dynldf/ ln_dynldf_lap  , ln_dynldf_bilap,                &
81         &                 ln_dynldf_level, ln_dynldf_hor, ln_dynldf_iso,   &
82         &                 ahm0, ahmb0
83      !!----------------------------------------------------------------------
84
85
86      ! Define the lateral physics parameters
87      ! ======================================
88   
89      ! Read Namelist nam_dynldf : Lateral physics
90      REWIND( numnam )
91      READ  ( numnam, nam_dynldf )
92
93      ! Parameter print
94      IF(lwp) THEN
95         WRITE(numout,*)
96         WRITE(numout,*) 'ldf_dyn : lateral momentum physics'
97         WRITE(numout,*) '~~~~~~~'
98         WRITE(numout,*) '          Namelist nam_dynldf : set lateral mixing parameters'
99         WRITE(numout,*) '             laplacian operator          ln_dynldf_lap   = ', ln_dynldf_lap
100         WRITE(numout,*) '             bilaplacian operator        ln_dynldf_bilap = ', ln_dynldf_bilap
101         WRITE(numout,*) '             iso-level                   ln_dynldf_level = ', ln_dynldf_level
102         WRITE(numout,*) '             horizontal (geopotential)   ln_dynldf_hor   = ', ln_dynldf_hor
103         WRITE(numout,*) '             iso-neutral                 ln_dynldf_iso   = ', ln_dynldf_iso
104         WRITE(numout,*) '             horizontal eddy viscosity            ahm0   = ', ahm0
105         WRITE(numout,*) '             background viscosity                 ahmb0  = ', ahmb0
106      ENDIF
107
108      ! ... check of lateral diffusive operator on tracers
109      !           ==> will be done in trazdf module
110
111      ! ... Space variation of eddy coefficients
112      ioptio = 0
113#if defined key_dynldf_c3d
114      IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude, depth)'
115      ioptio = ioptio+1
116#endif
117#if defined key_dynldf_c2d
118      IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude)'
119      ioptio = ioptio+1
120#endif
121#if defined key_dynldf_c1d
122      IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( depth )'
123      ioptio = ioptio+1
124      IF( ln_sco ) THEN
125         IF(lwp) WRITE(numout,cform_err)
126         IF(lwp) WRITE(numout,*) '          key_dynldf_c1d cannot be used in s-coordinate (ln_sco)'
127         nstop = nstop + 1
128      ENDIF
129#endif
130      IF( ioptio == 0 ) THEN
131          IF(lwp) WRITE(numout,*) '          momentum mixing coef. = constant  (default option)'
132        ELSEIF( ioptio > 1 ) THEN
133          IF(lwp) WRITE(numout,cform_err)
134          IF(lwp) WRITE(numout,*) '          use only one of the following keys:',   &
135                                  ' key_dynldf_c3d, key_dynldf_c2d, key_dynldf_c1d'
136          nstop = nstop + 1
137      ENDIF
138
139
140      IF( ln_dynldf_bilap ) THEN
141         IF(lwp) WRITE(numout,*) '          biharmonic momentum diffusion'
142         IF( ahm0 > 0 .AND. .NOT. lk_esopa ) THEN
143            IF(lwp) WRITE(numout,cform_err)
144            IF(lwp) WRITE(numout,*) 'The horizontal viscosity coef. ahm0 must be negative'
145            nstop = nstop + 1
146         ENDIF
147      ELSE
148         IF(lwp) WRITE(numout,*) '          harmonic momentum diff. (default)'
149         IF( ahm0 < 0 .AND. .NOT. lk_esopa ) THEN
150            IF(lwp) WRITE(numout,cform_err)
151            IF(lwp) WRITE(numout,*) '          The horizontal viscosity coef. ahm0 must be positive'
152            nstop = nstop + 1
153         ENDIF
154      ENDIF
155
156
157      ! Lateral eddy viscosity
158      ! ======================
159
160#if defined key_dynldf_c3d
161      CALL ldf_dyn_c3d( ll_print )   ! ahm = 3D coef. = F( longitude, latitude, depth )
162#elif defined key_dynldf_c2d
163      CALL ldf_dyn_c2d( ll_print )   ! ahm = 1D coef. = F( longitude, latitude )
164#elif defined key_dynldf_c1d
165      CALL ldf_dyn_c1d( ll_print )   ! ahm = 1D coef. = F( depth )
166#else
167                                     ! Constant coefficients
168      IF(lwp) WRITE(numout,*)
169      IF(lwp) WRITE(numout,*) 'inildf: constant eddy viscosity coef. '
170      IF(lwp) WRITE(numout,*) '~~~~~~'
171      IF(lwp) WRITE(numout,*) '        ahm1 = ahm2 = ahm0 =  ',ahm0
172#endif
173
174   END SUBROUTINE ldf_dyn_init
175
176#if defined key_dynldf_c3d
177#   include "ldfdyn_c3d.h90"
178#elif defined key_dynldf_c2d
179#   include "ldfdyn_c2d.h90"
180#elif defined key_dynldf_c1d
181#   include "ldfdyn_c1d.h90"
182#endif
183
184
185   SUBROUTINE ldf_zpf_1d( ld_print, pdam, pwam, pbot, pdep, pah )
186      !!----------------------------------------------------------------------
187      !!                  ***  ROUTINE ldf_zpf  ***
188      !!                   
189      !! ** Purpose :   vertical adimensional profile for eddy coefficient
190      !!
191      !! ** Method  :   1D eddy viscosity coefficients ( depth )
192      !!
193      !!----------------------------------------------------------------------
194      !! * Arguments
195      LOGICAL , INTENT (in   ) :: ld_print   ! If true, output arrays on numout
196      REAL(wp), INTENT (in   ) ::   &
197          pdam,     &  ! depth of the inflection point
198          pwam,     &  ! width of inflection
199          pbot         ! battom value (0<pbot<= 1)
200      REAL(wp), INTENT (in   ), DIMENSION(jpk) ::   &
201          pdep         ! depth of the gridpoint (T, U, V, F)
202      REAL(wp), INTENT (inout), DIMENSION(jpk) ::   &
203          pah          ! adimensional vertical profile
204
205      !! * Local variables
206      INTEGER  ::   jk           ! dummy loop indices
207      REAL(wp) ::   zm00, zm01, zmhb, zmhs       ! temporary scalars
208      !!----------------------------------------------------------------------
209
210      zm00 = TANH( ( pdam - gdept_0(1    ) ) / pwam )
211      zm01 = TANH( ( pdam - gdept_0(jpkm1) ) / pwam )
212      zmhs = zm00 / zm01
213      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
214
215      DO jk = 1, jpk
216         pah(jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam )  )
217      END DO
218
219      ! Control print
220      IF(lwp .AND. ld_print ) THEN
221         WRITE(numout,*)
222         WRITE(numout,*) '         ahm profile : '
223         WRITE(numout,*)
224         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
225         DO jk = 1, jpk
226            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(jk), pdep(jk)
227         END DO
228      ENDIF
229
230   END SUBROUTINE ldf_zpf_1d
231
232
233   SUBROUTINE ldf_zpf_1d_3d( ld_print, pdam, pwam, pbot, pdep, pah )
234      !!----------------------------------------------------------------------
235      !!                  ***  ROUTINE ldf_zpf  ***
236      !!
237      !! ** Purpose :   vertical adimensional profile for eddy coefficient
238      !!
239      !! ** Method  :   1D eddy viscosity coefficients ( depth )
240      !!
241      !!----------------------------------------------------------------------
242      !! * Arguments
243      LOGICAL , INTENT (in   ) :: ld_print   ! If true, output arrays on numout
244      REAL(wp), INTENT (in   ) ::   &
245          pdam,     &  ! depth of the inflection point
246          pwam,     &  ! width of inflection
247          pbot         ! battom value (0<pbot<= 1)
248      REAL(wp), INTENT (in   ), DIMENSION(jpk) ::   &
249          pdep         ! depth of the gridpoint (T, U, V, F)
250      REAL(wp), INTENT (inout), DIMENSION(jpi,jpj,jpk) ::   &
251          pah          ! adimensional vertical profile
252
253      !! * Local variables
254      INTEGER  ::   jk           ! dummy loop indices
255      REAL(wp) ::   zm00, zm01, zmhb, zmhs, zcf  ! temporary scalars
256      !!----------------------------------------------------------------------
257
258      zm00 = TANH( ( pdam - gdept_0(1    ) ) / pwam )
259      zm01 = TANH( ( pdam - gdept_0(jpkm1) ) / pwam )
260      zmhs = zm00 / zm01
261      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
262
263      DO jk = 1, jpk
264         zcf = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam )  )
265         pah(:,:,jk) = zcf
266      END DO
267
268      ! Control print
269      IF(lwp .AND. ld_print ) THEN
270         WRITE(numout,*)
271         WRITE(numout,*) '         ahm profile : '
272         WRITE(numout,*)
273         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
274         DO jk = 1, jpk
275            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(jk)
276         END DO
277      ENDIF
278
279   END SUBROUTINE ldf_zpf_1d_3d
280
281
282   SUBROUTINE ldf_zpf_3d( ld_print, pdam, pwam, pbot, pdep, pah )
283      !!----------------------------------------------------------------------
284      !!                  ***  ROUTINE ldf_zpf  ***
285      !!
286      !! ** Purpose :   vertical adimensional profile for eddy coefficient
287      !!
288      !! ** Method  :   3D for partial step or s-coordinate
289      !!
290      !!----------------------------------------------------------------------
291      !! * Arguments
292      LOGICAL , INTENT (in   ) :: ld_print   ! If true, output arrays on numout
293      REAL(wp), INTENT (in   ) ::   &
294          pdam,     &  ! depth of the inflection point
295          pwam,     &  ! width of inflection
296          pbot         ! reduction factor (surface value / bottom value)
297      REAL(wp), INTENT (in   ), DIMENSION(jpi,jpj,jpk) ::   &
298          pdep         ! dep of the gridpoint (T, U, V, F)
299      REAL(wp), INTENT (inout), DIMENSION(jpi,jpj,jpk) ::   &
300          pah          ! adimensional vertical profile
301
302      !! * Local variables
303      INTEGER  ::   jk           ! dummy loop indices
304      REAL(wp) ::   zm00, zm01, zmhb, zmhs       ! temporary scalars
305      !!----------------------------------------------------------------------
306
307      zm00 = TANH( ( pdam - gdept_0(1    ) ) / pwam )   
308      zm01 = TANH( ( pdam - gdept_0(jpkm1) ) / pwam )
309      zmhs = zm00 / zm01
310      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
311
312      DO jk = 1, jpk
313         pah(:,:,jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(:,:,jk) ) / pwam )  )
314      END DO
315
316      ! Control print
317      IF(lwp .AND. ld_print ) THEN
318         WRITE(numout,*)
319         WRITE(numout,*) '         ahm profile : '
320         WRITE(numout,*)
321         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
322         DO jk = 1, jpk
323            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(1,1,jk)
324         END DO
325      ENDIF
326
327   END SUBROUTINE ldf_zpf_3d
328   !!======================================================================
329END MODULE ldfdyn
Note: See TracBrowser for help on using the repository browser.