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 branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 5 years ago

The Dr Hook changes from my perl code.
File size: 15.3 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   USE yomhook, ONLY: lhook, dr_hook
28   USE parkind1, ONLY: jprb, jpim
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   ldf_dyn_init   ! called by opa.F90
34
35  INTERFACE ldf_zpf
36     MODULE PROCEDURE ldf_zpf_1d, ldf_zpf_1d_3d, ldf_zpf_3d
37  END INTERFACE
38
39   !! * Substitutions
40#  include "domzgr_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
43   !! $Id$
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE ldf_dyn_init
49      !!----------------------------------------------------------------------
50      !!                  ***  ROUTINE ldf_dyn_init  ***
51      !!                   
52      !! ** Purpose :   set the horizontal ocean dynamics physics
53      !!
54      !! ** Method  : 
55      !!      -  default option : ahm = constant coef. = rn_ahm_0 (namelist)
56      !!      - 'key_dynldf_c1d': ahm = F(depth)                     see ldf_dyn_c1d.h90
57      !!      - 'key_dynldf_c2d': ahm = F(latitude,longitude)        see ldf_dyn_c2d.h90
58      !!      - 'key_dynldf_c3d': ahm = F(latitude,longitude,depth)  see ldf_dyn_c3d.h90
59      !!
60      !!      N.B. User defined include files.  By default, 3d and 2d coef.
61      !!      are set to a constant value given in the namelist and the 1d
62      !!      coefficients are initialized to a hyperbolic tangent vertical
63      !!      profile.
64      !!
65      !! Reference :   Madec, G. and M. Imbard, 1996: Climate Dynamics, 12, 381-388.
66      !!----------------------------------------------------------------------
67      INTEGER ::   ioptio         ! ???
68      INTEGER ::   ios            ! Local : output status for namelist read
69      LOGICAL ::   ll_print = .FALSE.    ! Logical flag for printing viscosity coef.
70      !!
71      NAMELIST/namdyn_ldf/ ln_dynldf_lap  , ln_dynldf_bilap,                  &
72         &                 ln_dynldf_level, ln_dynldf_hor  , ln_dynldf_iso,   &
73         &                 rn_ahm_0_lap   , rn_ahmb_0      , rn_ahm_0_blp ,   &
74         &                 rn_cmsmag_1    , rn_cmsmag_2    , rn_cmsh,         &
75         &                 rn_ahm_m_lap   , rn_ahm_m_blp
76         INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
77         INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
78         REAL(KIND=jprb)               :: zhook_handle
79
80         CHARACTER(LEN=*), PARAMETER :: RoutineName='LDF_DYN_INIT'
81
82         IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
83
84
85   !!----------------------------------------------------------------------
86
87      REWIND( numnam_ref )              ! Namelist namdyn_ldf in reference namelist : Lateral physics
88      READ  ( numnam_ref, namdyn_ldf, IOSTAT = ios, ERR = 901)
89901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_ldf in reference namelist', lwp )
90
91      REWIND( numnam_cfg )              ! Namelist namdyn_ldf in configuration namelist : Lateral physics
92      READ  ( numnam_cfg, namdyn_ldf, IOSTAT = ios, ERR = 902 )
93902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_ldf in configuration namelist', lwp )
94      IF(lwm) WRITE ( numond, namdyn_ldf )
95
96      IF(lwp) THEN                      ! Parameter print
97         WRITE(numout,*)
98         WRITE(numout,*) 'ldf_dyn : lateral momentum physics'
99         WRITE(numout,*) '~~~~~~~'
100         WRITE(numout,*) '   Namelist namdyn_ldf : set lateral mixing parameters'
101         WRITE(numout,*) '      laplacian operator                      ln_dynldf_lap   = ', ln_dynldf_lap
102         WRITE(numout,*) '      bilaplacian operator                    ln_dynldf_bilap = ', ln_dynldf_bilap
103         WRITE(numout,*) '      iso-level                               ln_dynldf_level = ', ln_dynldf_level
104         WRITE(numout,*) '      horizontal (geopotential)               ln_dynldf_hor   = ', ln_dynldf_hor
105         WRITE(numout,*) '      iso-neutral                             ln_dynldf_iso   = ', ln_dynldf_iso
106         WRITE(numout,*) '      horizontal laplacian eddy viscosity     rn_ahm_0_lap    = ', rn_ahm_0_lap
107         WRITE(numout,*) '      background viscosity                    rn_ahmb_0       = ', rn_ahmb_0
108         WRITE(numout,*) '      horizontal bilaplacian eddy viscosity   rn_ahm_0_blp    = ', rn_ahm_0_blp
109         WRITE(numout,*) '      upper limit for laplacian eddy visc     rn_ahm_m_lap    = ', rn_ahm_m_lap
110         WRITE(numout,*) '      upper limit for bilap eddy viscosity    rn_ahm_m_blp    = ', rn_ahm_m_blp
111
112      ENDIF
113
114      ahm0     = rn_ahm_0_lap              ! OLD namelist variables defined from DOCTOR namelist variables
115      ahmb0    = rn_ahmb_0
116      ahm0_blp = rn_ahm_0_blp
117
118      ! ... check of lateral diffusive operator on tracers
119      !           ==> will be done in trazdf module
120
121      ! ... Space variation of eddy coefficients
122      ioptio = 0
123#if defined key_dynldf_c3d
124      IF(lwp) WRITE(numout,*) '   momentum mixing coef. = F( latitude, longitude, depth)'
125      ioptio = ioptio+1
126#endif
127#if defined key_dynldf_c2d
128      IF(lwp) WRITE(numout,*) '   momentum mixing coef. = F( latitude, longitude)'
129      ioptio = ioptio+1
130#endif
131#if defined key_dynldf_c1d
132      IF(lwp) WRITE(numout,*) '   momentum mixing coef. = F( depth )'
133      ioptio = ioptio+1
134      IF( ln_sco ) CALL ctl_stop( 'key_dynldf_c1d cannot be used in s-coordinate (ln_sco)' )
135#endif
136      IF( ioptio == 0 ) THEN
137          IF(lwp) WRITE(numout,*) '   momentum mixing coef. = constant  (default option)'
138        ELSEIF( ioptio > 1 ) THEN
139           CALL ctl_stop( 'use only one of the following keys: key_dynldf_c3d, key_dynldf_c2d, key_dynldf_c1d' )
140      ENDIF
141
142
143      IF( ln_dynldf_bilap ) THEN
144         IF(lwp) WRITE(numout,*) '   biharmonic momentum diffusion'
145         IF( .NOT. ln_dynldf_lap ) ahm0 = ahm0_blp   ! Allow spatially varying coefs, which use ahm0 as input
146         IF( ahm0_blp > 0 .AND. .NOT. lk_esopa )   CALL ctl_stop( 'The horizontal viscosity coef. ahm0 must be negative' )
147      ELSE
148         IF(lwp) WRITE(numout,*) '   harmonic momentum diff. (default)'
149         IF( ahm0 < 0 .AND. .NOT. lk_esopa )   CALL ctl_stop( 'The horizontal viscosity coef. ahm0 must be positive' )
150      ENDIF
151
152
153      ! Lateral eddy viscosity
154      ! ======================
155#if defined key_dynldf_c3d
156      CALL ldf_dyn_c3d( ll_print )   ! ahm = 3D coef. = F( longitude, latitude, depth )
157#elif defined key_dynldf_c2d
158      CALL ldf_dyn_c2d( ll_print )   ! ahm = 1D coef. = F( longitude, latitude )
159#elif defined key_dynldf_c1d
160      CALL ldf_dyn_c1d( ll_print )   ! ahm = 1D coef. = F( depth )
161#else
162                                     ! Constant coefficients
163      IF(lwp) WRITE(numout,*)
164      IF(lwp) WRITE(numout,*) 'inildf: constant eddy viscosity coef. '
165      IF(lwp) WRITE(numout,*) '~~~~~~'
166      IF(lwp) WRITE(numout,*) '        ahm1 = ahm2 = ahm0 =  ',ahm0
167#endif
168     nkahm_smag = 0
169#if defined key_dynldf_smag
170     nkahm_smag = 1
171#endif
172
173      !
174         IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
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      LOGICAL , INTENT(in   )                 ::   ld_print   ! If true, output arrays on numout
195      REAL(wp), INTENT(in   )                 ::   pdam       ! depth of the inflection point
196      REAL(wp), INTENT(in   )                 ::   pwam       ! width of inflection
197      REAL(wp), INTENT(in   )                 ::   pbot       ! bottom value (0<pbot<= 1)
198      REAL(wp), INTENT(in   ), DIMENSION(jpk) ::   pdep       ! depth of the gridpoint (T, U, V, F)
199      REAL(wp), INTENT(inout), DIMENSION(jpk) ::   pah        ! adimensional vertical profile
200      !!
201      INTEGER  ::   jk           ! dummy loop indices
202      REAL(wp) ::   zm00, zm01, zmhb, zmhs       ! temporary scalars
203      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
204      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
205      REAL(KIND=jprb)               :: zhook_handle
206
207      CHARACTER(LEN=*), PARAMETER :: RoutineName='LDF_ZPF_1D'
208
209      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
210
211      !!----------------------------------------------------------------------
212
213      zm00 = TANH( ( pdam - gdept_1d(1    ) ) / pwam )
214      zm01 = TANH( ( pdam - gdept_1d(jpkm1) ) / pwam )
215      zmhs = zm00 / zm01
216      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
217
218      DO jk = 1, jpk
219         pah(jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam )  )
220      END DO
221
222      IF(lwp .AND. ld_print ) THEN      ! Control print
223         WRITE(numout,*)
224         WRITE(numout,*) '         ahm profile : '
225         WRITE(numout,*)
226         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
227         DO jk = 1, jpk
228            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(jk), pdep(jk)
229         END DO
230      ENDIF
231      !
232      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
233   END SUBROUTINE ldf_zpf_1d
234
235
236   SUBROUTINE ldf_zpf_1d_3d( ld_print, pdam, pwam, pbot, pdep, pah )
237      !!----------------------------------------------------------------------
238      !!                  ***  ROUTINE ldf_zpf  ***
239      !!
240      !! ** Purpose :   vertical adimensional profile for eddy coefficient
241      !!
242      !! ** Method  :   1D eddy viscosity coefficients ( depth )
243      !!----------------------------------------------------------------------
244      LOGICAL , INTENT(in   )                         ::   ld_print   ! If true, output arrays on numout
245      REAL(wp), INTENT(in   )                         ::   pdam       ! depth of the inflection point
246      REAL(wp), INTENT(in   )                         ::   pwam       ! width of inflection
247      REAL(wp), INTENT(in   )                         ::   pbot       ! bottom value (0<pbot<= 1)
248      REAL(wp), INTENT(in   ), DIMENSION          (:) ::   pdep       ! depth of the gridpoint (T, U, V, F)
249      REAL(wp), INTENT(inout), DIMENSION      (:,:,:) ::   pah        ! adimensional vertical profile
250      !!
251      INTEGER  ::   jk           ! dummy loop indices
252      REAL(wp) ::   zm00, zm01, zmhb, zmhs, zcf  ! temporary scalars
253      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
254      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
255      REAL(KIND=jprb)               :: zhook_handle
256
257      CHARACTER(LEN=*), PARAMETER :: RoutineName='LDF_ZPF_1D_3D'
258
259      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
260
261      !!----------------------------------------------------------------------
262
263      zm00 = TANH( ( pdam - gdept_1d(1    ) ) / pwam )
264      zm01 = TANH( ( pdam - gdept_1d(jpkm1) ) / pwam )
265      zmhs = zm00 / zm01
266      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
267
268      DO jk = 1, jpk
269         zcf = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam )  )
270         pah(:,:,jk) = zcf
271      END DO
272
273      IF(lwp .AND. ld_print ) THEN      ! Control print
274         WRITE(numout,*)
275         WRITE(numout,*) '         ahm profile : '
276         WRITE(numout,*)
277         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
278         DO jk = 1, jpk
279            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(jk)
280         END DO
281      ENDIF
282      !
283      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
284   END SUBROUTINE ldf_zpf_1d_3d
285
286
287   SUBROUTINE ldf_zpf_3d( ld_print, pdam, pwam, pbot, pdep, pah )
288      !!----------------------------------------------------------------------
289      !!                  ***  ROUTINE ldf_zpf  ***
290      !!
291      !! ** Purpose :   vertical adimensional profile for eddy coefficient
292      !!
293      !! ** Method  :   3D for partial step or s-coordinate
294      !!----------------------------------------------------------------------
295      LOGICAL , INTENT(in   )                         ::   ld_print   ! If true, output arrays on numout
296      REAL(wp), INTENT(in   )                         ::   pdam       ! depth of the inflection point
297      REAL(wp), INTENT(in   )                         ::   pwam       ! width of inflection
298      REAL(wp), INTENT(in   )                         ::   pbot       ! bottom value (0<pbot<= 1)
299      REAL(wp), INTENT(in   ), DIMENSION      (:,:,:) ::   pdep       ! dep of the gridpoint (T, U, V, F)
300      REAL(wp), INTENT(inout), DIMENSION      (:,:,:) ::   pah        ! adimensional vertical profile
301      !!
302      INTEGER  ::   jk           ! dummy loop indices
303      REAL(wp) ::   zm00, zm01, zmhb, zmhs       ! temporary scalars
304      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
305      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
306      REAL(KIND=jprb)               :: zhook_handle
307
308      CHARACTER(LEN=*), PARAMETER :: RoutineName='LDF_ZPF_3D'
309
310      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
311
312      !!----------------------------------------------------------------------
313
314      zm00 = TANH( ( pdam - gdept_1d(1    ) ) / pwam )   
315      zm01 = TANH( ( pdam - gdept_1d(jpkm1) ) / pwam )
316      zmhs = zm00 / zm01
317      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
318
319      DO jk = 1, jpk
320         pah(:,:,jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(:,:,jk) ) / pwam )  )
321      END DO
322
323      IF(lwp .AND. ld_print ) THEN      ! Control print
324         WRITE(numout,*)
325         WRITE(numout,*) '         ahm profile : '
326         WRITE(numout,*)
327         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
328         DO jk = 1, jpk
329            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(1,1,jk)
330         END DO
331      ENDIF
332      !
333      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
334   END SUBROUTINE ldf_zpf_3d
335
336   !!======================================================================
337END MODULE ldfdyn
Note: See TracBrowser for help on using the repository browser.