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 on Ticket #371 – Attachment – NEMO

Ticket #371: ldfdyn.F90

File ldfdyn.F90, 13.6 KB (added by ed.blockley, 15 years ago)

LDF/ldfdyn.F90

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   !! $Id: ldfdyn.F90 1152 2008-06-26 14:11:13Z rblod $
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      !!      Bilaplacian and laplacian schemes can be used together.
64      !!      If both are used, bilap uses constant coefficient, ahm0_blp.
65      !!      If only bilap is used, its coefficient can vary spatially and
66      !!      if new variable ahm0_blp is not set by namelist, then ahm0 is
67      !!      used as before for backward compatability.
68      !!
69      !! Reference :
70      !!      Madec, G. and M. Imbard, 1996, A global ocean mesh to overcome
71      !!      the North Pole singularity, Climate Dynamics, 12, 381-388.
72      !!
73      !! History :
74      !!        !  07-97  (G. Madec)  from inimix.F split in 2 routines
75      !!        !  08-97  (G. Madec)  multi dimensional coefficients
76      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
77      !!----------------------------------------------------------------------
78      !! * Modules used
79      USE ioipsl
80
81      !! * Local declarations
82      INTEGER ::   ioptio         ! ???
83      LOGICAL :: ll_print = .FALSE.    ! Logical flag for printing viscosity coef.
84
85       
86      NAMELIST/nam_dynldf/ ln_dynldf_lap  , ln_dynldf_bilap,                &
87         &                 ln_dynldf_level, ln_dynldf_hor, ln_dynldf_iso,   &
88         &                 ahm0, ahmb0, ahm0_blp
89      !!----------------------------------------------------------------------
90
91
92      ! Define the lateral physics parameters
93      ! ======================================
94   
95      ! Read Namelist nam_dynldf : Lateral physics
96      REWIND( numnam )
97      READ  ( numnam, nam_dynldf )
98
99      ! Parameter print
100      IF(lwp) THEN
101         WRITE(numout,*)
102         WRITE(numout,*) 'ldf_dyn : lateral momentum physics'
103         WRITE(numout,*) '~~~~~~~'
104         WRITE(numout,*) '          Namelist nam_dynldf : set lateral mixing parameters'
105         WRITE(numout,*) '             laplacian operator          ln_dynldf_lap   = ', ln_dynldf_lap
106         WRITE(numout,*) '             bilaplacian operator        ln_dynldf_bilap = ', ln_dynldf_bilap
107         WRITE(numout,*) '             iso-level                   ln_dynldf_level = ', ln_dynldf_level
108         WRITE(numout,*) '             horizontal (geopotential)   ln_dynldf_hor   = ', ln_dynldf_hor
109         WRITE(numout,*) '             iso-neutral                 ln_dynldf_iso   = ', ln_dynldf_iso
110         WRITE(numout,*) '             horizontal laplacian eddy viscosity  ahm0   = ', ahm0
111         WRITE(numout,*) '             background laplacian viscosity       ahmb0  = ', ahmb0
112         WRITE(numout,*) '             horizontal eddy viscosity (bilap)  ahm0_blp = ', ahm0_blp
113      ENDIF
114
115      ! ... check of lateral diffusive operator on tracers
116      !           ==> will be done in trazdf module
117
118      ! ... Space variation of eddy coefficients
119      ioptio = 0
120#if defined key_dynldf_c3d
121      IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude, depth)'
122      ioptio = ioptio+1
123#endif
124#if defined key_dynldf_c2d
125      IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude)'
126      ioptio = ioptio+1
127#endif
128#if defined key_dynldf_c1d
129      IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( depth )'
130      ioptio = ioptio+1
131      IF( ln_sco ) CALL ctl_stop( '          key_dynldf_c1d cannot be used in s-coordinate (ln_sco)' )
132#endif
133      IF( ioptio == 0 ) THEN
134          IF(lwp) WRITE(numout,*) '          momentum mixing coef. = constant  (default option)'
135        ELSEIF( ioptio > 1 ) THEN
136           CALL ctl_stop( '          use only one of the following keys:',   &
137                &         ' key_dynldf_c3d, key_dynldf_c2d, key_dynldf_c1d' )
138      ENDIF
139
140
141      IF( ln_dynldf_bilap ) THEN
142         IF(lwp) WRITE(numout,*) '          biharmonic momentum diffusion'
143
144         IF( ahm0_blp == 0.0 ) ahm0_blp = ahm0      ! Old namelist method: bilap specified with ahm0
145
146         IF( .NOT. ln_dynldf_lap ) ahm0 = ahm0_blp   ! Allow spatially varying coefs, which use ahm0 as input
147
148         IF( ahm0_blp > 0 .AND. .NOT. lk_esopa )    &
149              &   CALL ctl_stop( 'The horizontal viscosity coef. ahm0 must be negative' )
150      ENDIF
151      IF( ln_dynldf_lap ) THEN
152         IF(lwp) WRITE(numout,*) '          harmonic momentum diff. (default)'
153         IF( ahm0 < 0 .AND. .NOT. lk_esopa ) &
154              &   CALL ctl_stop( '          The horizontal viscosity coef. ahm0 must be positive' )
155      ENDIF
156
157
158      ! Lateral eddy viscosity
159      ! ======================
160
161#if defined key_dynldf_c3d
162      CALL ldf_dyn_c3d( ll_print )   ! ahm = 3D coef. = F( longitude, latitude, depth )
163#elif defined key_dynldf_c2d
164      CALL ldf_dyn_c2d( ll_print )   ! ahm = 1D coef. = F( longitude, latitude )
165#elif defined key_dynldf_c1d
166      CALL ldf_dyn_c1d( ll_print )   ! ahm = 1D coef. = F( depth )
167#else
168                                     ! Constant coefficients
169      IF(lwp) WRITE(numout,*)
170      IF(lwp) WRITE(numout,*) 'inildf: constant eddy viscosity coef. '
171      IF(lwp) WRITE(numout,*) '~~~~~~'
172      IF(lwp) WRITE(numout,*) '        ahm1 = ahm2 = ahm0 =  ',ahm0
173#endif
174
175   END SUBROUTINE ldf_dyn_init
176
177#if defined key_dynldf_c3d
178#   include "ldfdyn_c3d.h90"
179#elif defined key_dynldf_c2d
180#   include "ldfdyn_c2d.h90"
181#elif defined key_dynldf_c1d
182#   include "ldfdyn_c1d.h90"
183#endif
184
185
186   SUBROUTINE ldf_zpf_1d( ld_print, pdam, pwam, pbot, pdep, pah )
187      !!----------------------------------------------------------------------
188      !!                  ***  ROUTINE ldf_zpf  ***
189      !!                   
190      !! ** Purpose :   vertical adimensional profile for eddy coefficient
191      !!
192      !! ** Method  :   1D eddy viscosity coefficients ( depth )
193      !!
194      !!----------------------------------------------------------------------
195      !! * Arguments
196      LOGICAL , INTENT (in   ) :: ld_print   ! If true, output arrays on numout
197      REAL(wp), INTENT (in   ) ::   &
198          pdam,     &  ! depth of the inflection point
199          pwam,     &  ! width of inflection
200          pbot         ! battom value (0<pbot<= 1)
201      REAL(wp), INTENT (in   ), DIMENSION(jpk) ::   &
202          pdep         ! depth of the gridpoint (T, U, V, F)
203      REAL(wp), INTENT (inout), DIMENSION(jpk) ::   &
204          pah          ! adimensional vertical profile
205
206      !! * Local variables
207      INTEGER  ::   jk           ! dummy loop indices
208      REAL(wp) ::   zm00, zm01, zmhb, zmhs       ! temporary scalars
209      !!----------------------------------------------------------------------
210
211      zm00 = TANH( ( pdam - gdept_0(1    ) ) / pwam )
212      zm01 = TANH( ( pdam - gdept_0(jpkm1) ) / pwam )
213      zmhs = zm00 / zm01
214      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
215
216      DO jk = 1, jpk
217         pah(jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam )  )
218      END DO
219
220      ! Control print
221      IF(lwp .AND. ld_print ) THEN
222         WRITE(numout,*)
223         WRITE(numout,*) '         ahm profile : '
224         WRITE(numout,*)
225         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
226         DO jk = 1, jpk
227            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(jk), pdep(jk)
228         END DO
229      ENDIF
230
231   END SUBROUTINE ldf_zpf_1d
232
233
234   SUBROUTINE ldf_zpf_1d_3d( ld_print, pdam, pwam, pbot, pdep, pah )
235      !!----------------------------------------------------------------------
236      !!                  ***  ROUTINE ldf_zpf  ***
237      !!
238      !! ** Purpose :   vertical adimensional profile for eddy coefficient
239      !!
240      !! ** Method  :   1D eddy viscosity coefficients ( depth )
241      !!
242      !!----------------------------------------------------------------------
243      !! * Arguments
244      LOGICAL , INTENT (in   ) :: ld_print   ! If true, output arrays on numout
245      REAL(wp), INTENT (in   ) ::   &
246          pdam,     &  ! depth of the inflection point
247          pwam,     &  ! width of inflection
248          pbot         ! battom value (0<pbot<= 1)
249      REAL(wp), INTENT (in   ), DIMENSION(jpk) ::   &
250          pdep         ! depth of the gridpoint (T, U, V, F)
251      REAL(wp), INTENT (inout), DIMENSION(jpi,jpj,jpk) ::   &
252          pah          ! adimensional vertical profile
253
254      !! * Local variables
255      INTEGER  ::   jk           ! dummy loop indices
256      REAL(wp) ::   zm00, zm01, zmhb, zmhs, zcf  ! temporary scalars
257      !!----------------------------------------------------------------------
258
259      zm00 = TANH( ( pdam - gdept_0(1    ) ) / pwam )
260      zm01 = TANH( ( pdam - gdept_0(jpkm1) ) / pwam )
261      zmhs = zm00 / zm01
262      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
263
264      DO jk = 1, jpk
265         zcf = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam )  )
266         pah(:,:,jk) = zcf
267      END DO
268
269      ! Control print
270      IF(lwp .AND. ld_print ) THEN
271         WRITE(numout,*)
272         WRITE(numout,*) '         ahm profile : '
273         WRITE(numout,*)
274         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
275         DO jk = 1, jpk
276            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(jk)
277         END DO
278      ENDIF
279
280   END SUBROUTINE ldf_zpf_1d_3d
281
282
283   SUBROUTINE ldf_zpf_3d( ld_print, pdam, pwam, pbot, pdep, pah )
284      !!----------------------------------------------------------------------
285      !!                  ***  ROUTINE ldf_zpf  ***
286      !!
287      !! ** Purpose :   vertical adimensional profile for eddy coefficient
288      !!
289      !! ** Method  :   3D for partial step or s-coordinate
290      !!
291      !!----------------------------------------------------------------------
292      !! * Arguments
293      LOGICAL , INTENT (in   ) :: ld_print   ! If true, output arrays on numout
294      REAL(wp), INTENT (in   ) ::   &
295          pdam,     &  ! depth of the inflection point
296          pwam,     &  ! width of inflection
297          pbot         ! reduction factor (surface value / bottom value)
298      REAL(wp), INTENT (in   ), DIMENSION(jpi,jpj,jpk) ::   &
299          pdep         ! dep of the gridpoint (T, U, V, F)
300      REAL(wp), INTENT (inout), DIMENSION(jpi,jpj,jpk) ::   &
301          pah          ! adimensional vertical profile
302
303      !! * Local variables
304      INTEGER  ::   jk           ! dummy loop indices
305      REAL(wp) ::   zm00, zm01, zmhb, zmhs       ! temporary scalars
306      !!----------------------------------------------------------------------
307
308      zm00 = TANH( ( pdam - gdept_0(1    ) ) / pwam )   
309      zm01 = TANH( ( pdam - gdept_0(jpkm1) ) / pwam )
310      zmhs = zm00 / zm01
311      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
312
313      DO jk = 1, jpk
314         pah(:,:,jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(:,:,jk) ) / pwam )  )
315      END DO
316
317      ! Control print
318      IF(lwp .AND. ld_print ) THEN
319         WRITE(numout,*)
320         WRITE(numout,*) '         ahm profile : '
321         WRITE(numout,*)
322         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
323         DO jk = 1, jpk
324            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(1,1,jk)
325         END DO
326      ENDIF
327
328   END SUBROUTINE ldf_zpf_3d
329   !!======================================================================
330END MODULE ldfdyn