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

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

#2050 first version. Compiled OK in moci test suite

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