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 @ 247

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

CL : Add CVS Header and CeCILL licence information
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 16.8 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         WRITE(numout,*)
107      ENDIF
108
109      ! Parameter control
110
111      ! control the input
112      ioptio = 0
113      IF( ln_dynldf_lap   )   ioptio = ioptio + 1
114      IF( ln_dynldf_bilap )   ioptio = ioptio + 1
115      IF( ioptio /= 1 )   THEN
116          IF(lwp) WRITE(numout,cform_err)
117          IF(lwp) WRITE(numout,*) '          use ONE of the 2 lap/bilap operator type on momentum'
118          nstop = nstop + 1
119      ENDIF
120      ioptio = 0
121      IF( ln_dynldf_level )   ioptio = ioptio + 1
122      IF( ln_dynldf_hor   )   ioptio = ioptio + 1
123      IF( ln_dynldf_iso   )   ioptio = ioptio + 1
124      IF( ioptio /= 1 ) THEN
125         IF(lwp) WRITE(numout,cform_err)
126         IF(lwp) WRITE(numout,*) '          use only ONE direction (level/hor/iso)'
127         nstop = nstop + 1
128      ENDIF
129
130      IF( lk_sco ) THEN          ! s-coordinates: rotation required for horizontal or isopycnal direction
131         IF( ( ln_dynldf_iso .OR. ln_dynldf_hor ) .AND. .NOT.lk_ldfslp ) THEN
132            IF(lwp) WRITE(numout,cform_err)
133            IF(lwp) WRITE(numout,*) '          the rotation of the viscous tensor require key_ldfslp'
134            IF( .NOT.lk_esopa )   nstop = nstop + 1
135         ENDIF
136      ELSE                       ! z-coordinates with/without partial step:
137         ln_dynldf_level = ln_dynldf_level .OR. ln_dynldf_hor      ! level mixing = horizontal mixing
138         ln_dynldf_hor   = .FALSE.
139         IF(lwp) WRITE(numout,*) '          horizontal mixing in z-coord or partial steps: force ln_dynldf_level = T'
140         IF(lwp) WRITE(numout,*) '                                                  and    force ln_dynldf_hor   = F'
141         IF( ln_dynldf_iso .AND. .NOT.lk_ldfslp ) THEN             ! rotation required for isopycnal mixing
142            IF(lwp) WRITE(numout,cform_err)
143            IF(lwp) WRITE(numout,*) '          the rotation of the viscous tensor require key_ldfslp'
144            IF( .NOT.lk_esopa )   nstop = nstop + 1
145         ENDIF
146      ENDIF
147
148      l_dynldf_lap     =       ln_dynldf_lap   .AND. ln_dynldf_level     ! iso-level   laplacian operator
149      l_dynldf_bilap   =       ln_dynldf_bilap .AND. ln_dynldf_level     ! iso-level bilaplacian operator
150      l_dynldf_bilapg  =       ln_dynldf_bilap .AND. ln_dynldf_hor       ! geopotential bilap. (s-coord)
151      l_dynldf_iso     =       ln_dynldf_lap   .AND.                  &  ! laplacian operator
152         &                   ( ln_dynldf_iso   .OR.  ln_dynldf_hor )     ! iso-neutral (z-coord) or horizontal (s-coord)
153
154      l_dynzdf_iso    = .FALSE.
155      IF( l_dynldf_iso )   l_dynzdf_iso = .TRUE.
156
157      ioptio = 0
158      IF( l_dynldf_lap     )   ioptio = ioptio + 1
159      IF( l_dynldf_bilap   )   ioptio = ioptio + 1
160      IF( l_dynldf_bilapg  )   ioptio = ioptio + 1
161      IF( l_dynldf_iso     )   ioptio = ioptio + 1
162      IF( ioptio /= 1 ) THEN
163         IF(lwp) WRITE(numout,cform_err)
164         IF(lwp) WRITE(numout,*) '          this combination of operator and direction has not been implemented'
165         nstop = nstop + 1
166      ENDIF
167      IF( lk_esopa ) THEN
168         l_dynldf_lap = .TRUE.   ;   l_dynldf_bilap   = .TRUE.   ;   l_dynldf_bilapg  = .TRUE.
169         l_dynldf_iso = .TRUE.   ;   l_dynzdf_iso     = .TRUE.
170         IF(lwp ) WRITE(numout,*) '          esopa test: use all lateral physics options'
171      ENDIF
172
173      ! control print
174      IF( l_dynldf_lap    .AND. lwp ) WRITE(numout,*) '          iso-level laplacian momentum operator'
175      IF( l_dynldf_bilap  .AND. lwp ) WRITE(numout,*) '          iso-level bilaplacian momentum operator'
176      IF( l_dynldf_bilapg .AND. lwp ) WRITE(numout,*) '          geopotential bilaplacian momentum operator'
177      IF( l_dynldf_iso    .AND. lwp ) WRITE(numout,*) '          iso-neutral laplacian momentum operator'
178
179      ! ... Space variation of eddy coefficients
180      ioptio = 0
181#if defined key_dynldf_c3d
182      IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude, depth)'
183      ioptio = ioptio+1
184#endif
185#if defined key_dynldf_c2d
186      IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude)'
187      ioptio = ioptio+1
188#endif
189#if defined key_dynldf_c1d
190      IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( depth )'
191      ioptio = ioptio+1
192      IF( lk_sco ) THEN
193         IF(lwp) WRITE(numout,cform_err)
194         IF(lwp) WRITE(numout,*) '          key_dynldf_c1d cannot be used in s-coordinate (key_s_coord)'
195         nstop = nstop + 1
196      ENDIF
197#endif
198      IF( ioptio == 0 ) THEN
199          IF(lwp) WRITE(numout,*) '          momentum mixing coef. = constant  (default option)'
200        ELSEIF( ioptio > 1 ) THEN
201          IF(lwp) WRITE(numout,cform_err)
202          IF(lwp) WRITE(numout,*) '          use only one of the following keys:',   &
203                                  ' key_dynldf_c3d, key_dynldf_c2d, key_dynldf_c1d'
204          nstop = nstop + 1
205      ENDIF
206
207
208      IF( l_dynldf_bilap .OR. l_dynldf_bilapg ) THEN
209         IF(lwp) WRITE(numout,*) '          biharmonic momentum diffusion'
210         IF( ahm0 > 0 .AND. .NOT. lk_esopa ) THEN
211            IF(lwp) WRITE(numout,cform_err)
212            IF(lwp) WRITE(numout,*) 'The horizontal viscosity coef. ahm0 must be negative'
213            nstop = nstop + 1
214         ENDIF
215      ELSE
216         IF(lwp) WRITE(numout,*) '          harmonic momentum diff. (default)'
217         IF( ahm0 < 0 .AND. .NOT. lk_esopa ) THEN
218            IF(lwp) WRITE(numout,cform_err)
219            IF(lwp) WRITE(numout,*) '          The horizontal viscosity coef. ahm0 must be positive'
220            nstop = nstop + 1
221         ENDIF
222      ENDIF
223
224
225      ! Lateral eddy viscosity
226      ! ======================
227
228#if defined key_dynldf_c3d
229      CALL ldf_dyn_c3d( ll_print )   ! ahm = 3D coef. = F( longitude, latitude, depth )
230#elif defined key_dynldf_c2d
231      CALL ldf_dyn_c2d( ll_print )   ! ahm = 1D coef. = F( longitude, latitude )
232#elif defined key_dynldf_c1d
233      CALL ldf_dyn_c1d( ll_print )   ! ahm = 1D coef. = F( depth )
234#else
235                                     ! Constant coefficients
236      IF(lwp) WRITE(numout,*)
237      IF(lwp) WRITE(numout,*) 'inildf: constant eddy viscosity coef. '
238      IF(lwp) WRITE(numout,*) '~~~~~~'
239      IF(lwp) WRITE(numout,*) '        ahm1 = ahm2 = ahm0 =  ',ahm0
240#endif
241
242   END SUBROUTINE ldf_dyn_init
243
244#if defined key_dynldf_c3d
245#   include "ldfdyn_c3d.h90"
246#elif defined key_dynldf_c2d
247#   include "ldfdyn_c2d.h90"
248#elif defined key_dynldf_c1d
249#   include "ldfdyn_c1d.h90"
250#endif
251
252
253   SUBROUTINE ldf_zpf_1d( ld_print, pdam, pwam, pbot, pdep, pah )
254      !!----------------------------------------------------------------------
255      !!                  ***  ROUTINE ldf_zpf  ***
256      !!                   
257      !! ** Purpose :   vertical adimensional profile for eddy coefficient
258      !!
259      !! ** Method  :   1D eddy viscosity coefficients ( depth )
260      !!
261      !!----------------------------------------------------------------------
262      !! * Arguments
263      LOGICAL , INTENT (in   ) :: ld_print   ! If true, output arrays on numout
264      REAL(wp), INTENT (in   ) ::   &
265          pdam,     &  ! depth of the inflection point
266          pwam,     &  ! width of inflection
267          pbot         ! battom value (0<pbot<= 1)
268      REAL(wp), INTENT (in   ), DIMENSION(jpk) ::   &
269          pdep         ! depth of the gridpoint (T, U, V, F)
270      REAL(wp), INTENT (inout), DIMENSION(jpk) ::   &
271          pah          ! adimensional vertical profile
272
273      !! * Local variables
274      INTEGER  ::   jk           ! dummy loop indices
275      REAL(wp) ::   zm00, zm01, zmhb, zmhs       ! temporary scalars
276      !!----------------------------------------------------------------------
277
278      zm00 = TANH( ( pdam - gdept(1    ) ) / pwam )
279      zm01 = TANH( ( pdam - gdept(jpkm1) ) / pwam )
280      zmhs = zm00 / zm01
281      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
282
283      DO jk = 1, jpk
284         pah(jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam )  )
285      END DO
286
287      ! Control print
288      IF(lwp .AND. ld_print ) THEN
289         WRITE(numout,*)
290         WRITE(numout,*) '         ahm profile : '
291         WRITE(numout,*)
292         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
293         DO jk = 1, jpk
294            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(jk), pdep(jk)
295         END DO
296      ENDIF
297
298   END SUBROUTINE ldf_zpf_1d
299
300
301   SUBROUTINE ldf_zpf_1d_3d( ld_print, pdam, pwam, pbot, pdep, pah )
302      !!----------------------------------------------------------------------
303      !!                  ***  ROUTINE ldf_zpf  ***
304      !!
305      !! ** Purpose :   vertical adimensional profile for eddy coefficient
306      !!
307      !! ** Method  :   1D eddy viscosity coefficients ( depth )
308      !!
309      !!----------------------------------------------------------------------
310      !! * Arguments
311      LOGICAL , INTENT (in   ) :: ld_print   ! If true, output arrays on numout
312      REAL(wp), INTENT (in   ) ::   &
313          pdam,     &  ! depth of the inflection point
314          pwam,     &  ! width of inflection
315          pbot         ! battom value (0<pbot<= 1)
316      REAL(wp), INTENT (in   ), DIMENSION(jpk) ::   &
317          pdep         ! depth of the gridpoint (T, U, V, F)
318      REAL(wp), INTENT (inout), DIMENSION(jpi,jpj,jpk) ::   &
319          pah          ! adimensional vertical profile
320
321      !! * Local variables
322      INTEGER  ::   jk           ! dummy loop indices
323      REAL(wp) ::   zm00, zm01, zmhb, zmhs, zcf  ! temporary scalars
324      !!----------------------------------------------------------------------
325
326      zm00 = TANH( ( pdam - gdept(1    ) ) / pwam )
327      zm01 = TANH( ( pdam - gdept(jpkm1) ) / pwam )
328      zmhs = zm00 / zm01
329      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
330
331      DO jk = 1, jpk
332         zcf = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam )  )
333         pah(:,:,jk) = zcf
334      END DO
335
336      ! Control print
337      IF(lwp .AND. ld_print ) THEN
338         WRITE(numout,*)
339         WRITE(numout,*) '         ahm profile : '
340         WRITE(numout,*)
341         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
342         DO jk = 1, jpk
343            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(jk)
344         END DO
345      ENDIF
346
347   END SUBROUTINE ldf_zpf_1d_3d
348
349
350   SUBROUTINE ldf_zpf_3d( ld_print, pdam, pwam, pbot, pdep, pah )
351      !!----------------------------------------------------------------------
352      !!                  ***  ROUTINE ldf_zpf  ***
353      !!
354      !! ** Purpose :   vertical adimensional profile for eddy coefficient
355      !!
356      !! ** Method  :   3D for partial step or s-coordinate
357      !!
358      !!----------------------------------------------------------------------
359      !! * Arguments
360      LOGICAL , INTENT (in   ) :: ld_print   ! If true, output arrays on numout
361      REAL(wp), INTENT (in   ) ::   &
362          pdam,     &  ! depth of the inflection point
363          pwam,     &  ! width of inflection
364          pbot         ! reduction factor (surface value / bottom value)
365      REAL(wp), INTENT (in   ), DIMENSION(jpi,jpj,jpk) ::   &
366          pdep         ! dep of the gridpoint (T, U, V, F)
367      REAL(wp), INTENT (inout), DIMENSION(jpi,jpj,jpk) ::   &
368          pah          ! adimensional vertical profile
369
370      !! * Local variables
371      INTEGER  ::   jk           ! dummy loop indices
372      REAL(wp) ::   zm00, zm01, zmhb, zmhs       ! temporary scalars
373      !!----------------------------------------------------------------------
374
375      zm00 = TANH( ( pdam - gdept(1    ) ) / pwam )   
376      zm01 = TANH( ( pdam - gdept(jpkm1) ) / pwam )
377      zmhs = zm00 / zm01
378      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
379
380      DO jk = 1, jpk
381         pah(:,:,jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(:,:,jk) ) / pwam )  )
382      END DO
383
384      ! Control print
385      IF(lwp .AND. ld_print ) THEN
386         WRITE(numout,*)
387         WRITE(numout,*) '         ahm profile : '
388         WRITE(numout,*)
389         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
390         DO jk = 1, jpk
391            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(1,1,jk)
392         END DO
393      ENDIF
394
395   END SUBROUTINE ldf_zpf_3d
396   !!======================================================================
397END MODULE ldfdyn
Note: See TracBrowser for help on using the repository browser.