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

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 9383

Last change on this file since 9383 was 9383, checked in by andmirek, 6 years ago

#2050 fixes and changes

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