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

source: tags/nemo_dev_x3/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 105

Last change on this file since 105 was 105, checked in by cvs2svn, 20 years ago

This commit was manufactured by cvs2svn to create tag 'nemo_dev_x3'.
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 16.7 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 , LODYC-IPSL  (2003)
37   !!----------------------------------------------------------------------
38
39CONTAINS
40
41   SUBROUTINE ldf_dyn_init
42      !!----------------------------------------------------------------------
43      !!                  ***  ROUTINE ldf_dyn_init  ***
44      !!                   
45      !! ** Purpose :   set the horizontal ocean dynamics physics
46      !!
47      !! ** Method  : 
48      !!      Eddy viscosity coefficients:
49      !!         default option   : constant coef. ahm0 (namelist)
50      !!        'key_dynldf_c1d': depth dependent coef. defined in
51      !!                        in ldf_dyn_c1d routine
52      !!        'key_dynldf_c2d': latitude and longitude dependent coef.
53      !!                        defined in ldf_dyn_c2d routine
54      !!        'key_dynldf_c3d': latitude, longitude, depth dependent coef.
55      !!                        defined in ldf_dyn_c3d routine
56      !!      N.B. User defined include files.  By default, 3d and 2d coef.
57      !!      are set to a constant value given in the namelist and the 1d
58      !!      coefficients are initialized to a hyperbolic tangent vertical
59      !!      profile.
60      !!
61      !! Reference :
62      !!      Madec, G. and M. Imbard, 1996, A global ocean mesh to overcome
63      !!      the North Pole singularity, Climate Dynamics, 12, 381-388.
64      !!
65      !! History :
66      !!        !  07-97  (G. Madec)  from inimix.F split in 2 routines
67      !!        !  08-97  (G. Madec)  multi dimensional coefficients
68      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
69      !!----------------------------------------------------------------------
70      !! * Modules used
71      USE ioipsl
72
73      !! * Local declarations
74      INTEGER ::   ioptio         ! ???
75      LOGICAL :: ll_print = .FALSE.    ! Logical flag for printing viscosity coef.
76
77       
78      NAMELIST/nam_dynldf/ ln_dynldf_lap  , ln_dynldf_bilap,                &
79         &                 ln_dynldf_level, ln_dynldf_hor, ln_dynldf_iso,   &
80         &                 ahm0, ahmb0
81      !!----------------------------------------------------------------------
82
83
84      ! Define the lateral physics parameters
85      ! ======================================
86   
87      ! Read Namelist nam_dynldf : Lateral physics
88      REWIND( numnam )
89      READ  ( numnam, nam_dynldf )
90
91      ! Parameter print
92      IF(lwp) THEN
93         WRITE(numout,*)
94         WRITE(numout,*) 'ldf_dyn : lateral momentum physics'
95         WRITE(numout,*) '~~~~~~~'
96         WRITE(numout,*) '          Namelist nam_dynldf : set lateral mixing parameters'
97         WRITE(numout,*) '             laplacian operator          ln_dynldf_lap   = ', ln_dynldf_lap
98         WRITE(numout,*) '             bilaplacian operator        ln_dynldf_bilap = ', ln_dynldf_bilap
99         WRITE(numout,*) '             iso-level                   ln_dynldf_level = ', ln_dynldf_level
100         WRITE(numout,*) '             horizontal (geopotential)   ln_dynldf_hor   = ', ln_dynldf_hor
101         WRITE(numout,*) '             iso-neutral                 ln_dynldf_iso   = ', ln_dynldf_iso
102         WRITE(numout,*) '             horizontal eddy viscosity            ahm0   = ', ahm0
103         WRITE(numout,*) '             background viscosity                 ahmb0  = ', ahmb0
104         WRITE(numout,*)
105      ENDIF
106
107      ! Parameter control
108
109      ! control the input
110      ioptio = 0
111      IF( ln_dynldf_lap   )   ioptio = ioptio + 1
112      IF( ln_dynldf_bilap )   ioptio = ioptio + 1
113      IF( ioptio /= 1 )   THEN
114          IF(lwp) WRITE(numout,cform_err)
115          IF(lwp) WRITE(numout,*) '          use ONE of the 2 lap/bilap operator type on momentum'
116          nstop = nstop + 1
117      ENDIF
118      ioptio = 0
119      IF( ln_dynldf_level )   ioptio = ioptio + 1
120      IF( ln_dynldf_hor   )   ioptio = ioptio + 1
121      IF( ln_dynldf_iso   )   ioptio = ioptio + 1
122      IF( ioptio /= 1 ) THEN
123         IF(lwp) WRITE(numout,cform_err)
124         IF(lwp) WRITE(numout,*) '          use only ONE direction (level/hor/iso)'
125         nstop = nstop + 1
126      ENDIF
127
128      IF( lk_sco ) THEN          ! s-coordinates: rotation required for horizontal or isopycnal direction
129         IF( ( ln_dynldf_iso .OR. ln_dynldf_hor ) .AND. .NOT.lk_ldfslp ) THEN
130            IF(lwp) WRITE(numout,cform_err)
131            IF(lwp) WRITE(numout,*) '          the rotation of the viscous tensor require key_ldfslp'
132            IF( .NOT.lk_esopa )   nstop = nstop + 1
133         ENDIF
134      ELSE                       ! z-coordinates with/without partial step:
135         ln_dynldf_level = ln_dynldf_level .OR. ln_dynldf_hor      ! level mixing = horizontal mixing
136         ln_dynldf_hor   = .FALSE.
137         IF(lwp) WRITE(numout,*) '          horizontal mixing in z-coord or partial steps: force ln_dynldf_level = T'
138         IF(lwp) WRITE(numout,*) '                                                  and    force ln_dynldf_hor   = F'
139         IF( ln_dynldf_iso .AND. .NOT.lk_ldfslp ) THEN             ! rotation required for isopycnal mixing
140            IF(lwp) WRITE(numout,cform_err)
141            IF(lwp) WRITE(numout,*) '          the rotation of the viscous tensor require key_ldfslp'
142            IF( .NOT.lk_esopa )   nstop = nstop + 1
143         ENDIF
144      ENDIF
145
146      l_dynldf_lap     =       ln_dynldf_lap   .AND. ln_dynldf_level     ! iso-level   laplacian operator
147      l_dynldf_bilap   =       ln_dynldf_bilap .AND. ln_dynldf_level     ! iso-level bilaplacian operator
148      l_dynldf_bilapg  =       ln_dynldf_bilap .AND. ln_dynldf_hor       ! geopotential bilap. (s-coord)
149      l_dynldf_iso     =       ln_dynldf_lap   .AND.                  &  ! laplacian operator
150         &                   ( ln_dynldf_iso   .OR.  ln_dynldf_hor )     ! iso-neutral (z-coord) or horizontal (s-coord)
151
152      l_dynzdf_iso    = .FALSE.
153      IF( l_dynldf_iso )   l_dynzdf_iso = .TRUE.
154
155      ioptio = 0
156      IF( l_dynldf_lap     )   ioptio = ioptio + 1
157      IF( l_dynldf_bilap   )   ioptio = ioptio + 1
158      IF( l_dynldf_bilapg  )   ioptio = ioptio + 1
159      IF( l_dynldf_iso     )   ioptio = ioptio + 1
160      IF( ioptio /= 1 ) THEN
161         IF(lwp) WRITE(numout,cform_err)
162         IF(lwp) WRITE(numout,*) '          this combination of operator and direction has not been implemented'
163         nstop = nstop + 1
164      ENDIF
165      IF( lk_esopa ) THEN
166         l_dynldf_lap = .TRUE.   ;   l_dynldf_bilap   = .TRUE.   ;   l_dynldf_bilapg  = .TRUE.
167         l_dynldf_iso = .TRUE.   ;   l_dynzdf_iso     = .TRUE.
168         IF(lwp ) WRITE(numout,*) '          esopa test: use all lateral physics options'
169      ENDIF
170
171      ! control print
172      IF( l_dynldf_lap    .AND. lwp ) WRITE(numout,*) '          iso-level laplacian momentum operator'
173      IF( l_dynldf_bilap  .AND. lwp ) WRITE(numout,*) '          iso-level bilaplacian momentum operator'
174      IF( l_dynldf_bilapg .AND. lwp ) WRITE(numout,*) '          geopotential bilaplacian momentum operator'
175      IF( l_dynldf_iso    .AND. lwp ) WRITE(numout,*) '          iso-neutral laplacian momentum operator'
176
177      ! ... Space variation of eddy coefficients
178      ioptio = 0
179#if defined key_dynldf_c3d
180      IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude, depth)'
181      ioptio = ioptio+1
182#endif
183#if defined key_dynldf_c2d
184      IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude)'
185      ioptio = ioptio+1
186#endif
187#if defined key_dynldf_c1d
188      IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( depth )'
189      ioptio = ioptio+1
190      IF( lk_sco ) THEN
191         IF(lwp) WRITE(numout,cform_err)
192         IF(lwp) WRITE(numout,*) '          key_dynldf_c1d cannot be used in s-coordinate (key_s_coord)'
193         nstop = nstop + 1
194      ENDIF
195#endif
196      IF( ioptio == 0 ) THEN
197          IF(lwp) WRITE(numout,*) '          momentum mixing coef. = constant  (default option)'
198        ELSEIF( ioptio > 1 ) THEN
199          IF(lwp) WRITE(numout,cform_err)
200          IF(lwp) WRITE(numout,*) '          use only one of the following keys:',   &
201                                  ' key_dynldf_c3d, key_dynldf_c2d, key_dynldf_c1d'
202          nstop = nstop + 1
203      ENDIF
204
205
206      IF( l_dynldf_bilap .OR. l_dynldf_bilapg ) THEN
207         IF(lwp) WRITE(numout,*) '          biharmonic momentum diffusion'
208         IF( ahm0 > 0 .AND. .NOT. lk_esopa ) THEN
209            IF(lwp) WRITE(numout,cform_err)
210            IF(lwp) WRITE(numout,*) 'The horizontal viscosity coef. ahm0 must be negative'
211            nstop = nstop + 1
212         ENDIF
213      ELSE
214         IF(lwp) WRITE(numout,*) '          harmonic momentum diff. (default)'
215         IF( ahm0 < 0 .AND. .NOT. lk_esopa ) THEN
216            IF(lwp) WRITE(numout,cform_err)
217            IF(lwp) WRITE(numout,*) '          The horizontal viscosity coef. ahm0 must be positive'
218            nstop = nstop + 1
219         ENDIF
220      ENDIF
221
222
223      ! Lateral eddy viscosity
224      ! ======================
225
226#if defined key_dynldf_c3d
227      CALL ldf_dyn_c3d( ll_print )   ! ahm = 3D coef. = F( longitude, latitude, depth )
228#elif defined key_dynldf_c2d
229      CALL ldf_dyn_c2d( ll_print )   ! ahm = 1D coef. = F( longitude, latitude )
230#elif defined key_dynldf_c1d
231      CALL ldf_dyn_c1d( ll_print )   ! ahm = 1D coef. = F( depth )
232#else
233                                     ! Constant coefficients
234      IF(lwp) WRITE(numout,*)
235      IF(lwp) WRITE(numout,*) 'inildf: constant eddy viscosity coef. '
236      IF(lwp) WRITE(numout,*) '~~~~~~'
237      IF(lwp) WRITE(numout,*) '        ahm1 = ahm2 = ahm0 =  ',ahm0
238#endif
239
240   END SUBROUTINE ldf_dyn_init
241
242#if defined key_dynldf_c3d
243#   include "ldfdyn_c3d.h90"
244#elif defined key_dynldf_c2d
245#   include "ldfdyn_c2d.h90"
246#elif defined key_dynldf_c1d
247#   include "ldfdyn_c1d.h90"
248#endif
249
250
251   SUBROUTINE ldf_zpf_1d( ld_print, pdam, pwam, pbot, pdep, pah )
252      !!----------------------------------------------------------------------
253      !!                  ***  ROUTINE ldf_zpf  ***
254      !!                   
255      !! ** Purpose :   vertical adimensional profile for eddy coefficient
256      !!
257      !! ** Method  :   1D eddy viscosity coefficients ( depth )
258      !!
259      !!----------------------------------------------------------------------
260      !! * Arguments
261      LOGICAL , INTENT (in   ) :: ld_print   ! If true, output arrays on numout
262      REAL(wp), INTENT (in   ) ::   &
263          pdam,     &  ! depth of the inflection point
264          pwam,     &  ! width of inflection
265          pbot         ! battom value (0<pbot<= 1)
266      REAL(wp), INTENT (in   ), DIMENSION(jpk) ::   &
267          pdep         ! depth of the gridpoint (T, U, V, F)
268      REAL(wp), INTENT (inout), DIMENSION(jpk) ::   &
269          pah          ! adimensional vertical profile
270
271      !! * Local variables
272      INTEGER  ::   jk           ! dummy loop indices
273      REAL(wp) ::   zm00, zm01, zmhb, zmhs       ! temporary scalars
274      !!----------------------------------------------------------------------
275
276      zm00 = TANH( ( pdam - gdept(1    ) ) / pwam )
277      zm01 = TANH( ( pdam - gdept(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      ! Control print
286      IF(lwp .AND. ld_print ) THEN
287         WRITE(numout,*)
288         WRITE(numout,*) '         ahm profile : '
289         WRITE(numout,*)
290         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
291         DO jk = 1, jpk
292            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(jk), pdep(jk)
293         END DO
294      ENDIF
295
296   END SUBROUTINE ldf_zpf_1d
297
298
299   SUBROUTINE ldf_zpf_1d_3d( ld_print, pdam, pwam, pbot, pdep, pah )
300      !!----------------------------------------------------------------------
301      !!                  ***  ROUTINE ldf_zpf  ***
302      !!
303      !! ** Purpose :   vertical adimensional profile for eddy coefficient
304      !!
305      !! ** Method  :   1D eddy viscosity coefficients ( depth )
306      !!
307      !!----------------------------------------------------------------------
308      !! * Arguments
309      LOGICAL , INTENT (in   ) :: ld_print   ! If true, output arrays on numout
310      REAL(wp), INTENT (in   ) ::   &
311          pdam,     &  ! depth of the inflection point
312          pwam,     &  ! width of inflection
313          pbot         ! battom value (0<pbot<= 1)
314      REAL(wp), INTENT (in   ), DIMENSION(jpk) ::   &
315          pdep         ! depth of the gridpoint (T, U, V, F)
316      REAL(wp), INTENT (inout), DIMENSION(jpi,jpj,jpk) ::   &
317          pah          ! adimensional vertical profile
318
319      !! * Local variables
320      INTEGER  ::   jk           ! dummy loop indices
321      REAL(wp) ::   zm00, zm01, zmhb, zmhs, zcf  ! temporary scalars
322      !!----------------------------------------------------------------------
323
324      zm00 = TANH( ( pdam - gdept(1    ) ) / pwam )
325      zm01 = TANH( ( pdam - gdept(jpkm1) ) / pwam )
326      zmhs = zm00 / zm01
327      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
328
329      DO jk = 1, jpk
330         zcf = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam )  )
331         pah(:,:,jk) = zcf
332      END DO
333
334      ! Control print
335      IF(lwp .AND. ld_print ) THEN
336         WRITE(numout,*)
337         WRITE(numout,*) '         ahm profile : '
338         WRITE(numout,*)
339         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
340         DO jk = 1, jpk
341            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(:,:,jk), pdep(jk)
342         END DO
343      ENDIF
344
345   END SUBROUTINE ldf_zpf_1d_3d
346
347
348   SUBROUTINE ldf_zpf_3d( ld_print, pdam, pwam, pbot, pdep, pah )
349      !!----------------------------------------------------------------------
350      !!                  ***  ROUTINE ldf_zpf  ***
351      !!
352      !! ** Purpose :   vertical adimensional profile for eddy coefficient
353      !!
354      !! ** Method  :   3D for partial step or s-coordinate
355      !!
356      !!----------------------------------------------------------------------
357      !! * Arguments
358      LOGICAL , INTENT (in   ) :: ld_print   ! If true, output arrays on numout
359      REAL(wp), INTENT (in   ) ::   &
360          pdam,     &  ! depth of the inflection point
361          pwam,     &  ! width of inflection
362          pbot         ! reduction factor (surface value / bottom value)
363      REAL(wp), INTENT (in   ), DIMENSION(jpi,jpj,jpk) ::   &
364          pdep         ! dep of the gridpoint (T, U, V, F)
365      REAL(wp), INTENT (inout), DIMENSION(jpi,jpj,jpk) ::   &
366          pah          ! adimensional vertical profile
367
368      !! * Local variables
369      INTEGER  ::   jk           ! dummy loop indices
370      REAL(wp) ::   zm00, zm01, zmhb, zmhs       ! temporary scalars
371      !!----------------------------------------------------------------------
372
373      zm00 = TANH( ( pdam - gdept(1    ) ) / pwam )   
374      zm01 = TANH( ( pdam - gdept(jpkm1) ) / pwam )
375      zmhs = zm00 / zm01
376      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
377
378      DO jk = 1, jpk
379         pah(:,:,jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(:,:,jk) ) / pwam )  )
380      END DO
381
382      ! Control print
383      IF(lwp .AND. ld_print ) THEN
384         WRITE(numout,*)
385         WRITE(numout,*) '         ahm profile : '
386         WRITE(numout,*)
387         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
388         DO jk = 1, jpk
389            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(1,1,jk)
390         END DO
391      ENDIF
392
393   END SUBROUTINE ldf_zpf_3d
394   !!======================================================================
395END MODULE ldfdyn
Note: See TracBrowser for help on using the repository browser.