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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 2690

Last change on this file since 2690 was 2690, checked in by gm, 13 years ago

dynamic mem: #785 ; homogeneization of the coding style associated with dyn allocation
  • Property svn:keywords set to Id
File size: 12.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
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   ldf_dyn_init   ! called by opa.F90
30
31  INTERFACE ldf_zpf
32     MODULE PROCEDURE ldf_zpf_1d, ldf_zpf_1d_3d, ldf_zpf_3d
33  END INTERFACE
34
35   !! * Substitutions
36#  include "domzgr_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
39   !! $Id$
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE ldf_dyn_init
45      !!----------------------------------------------------------------------
46      !!                  ***  ROUTINE ldf_dyn_init  ***
47      !!                   
48      !! ** Purpose :   set the horizontal ocean dynamics physics
49      !!
50      !! ** Method  : 
51      !!      -  default option : ahm = constant coef. = rn_ahm_0 (namelist)
52      !!      - 'key_dynldf_c1d': ahm = F(depth)                     see ldf_dyn_c1d.h90
53      !!      - 'key_dynldf_c2d': ahm = F(latitude,longitude)        see ldf_dyn_c2d.h90
54      !!      - 'key_dynldf_c3d': ahm = F(latitude,longitude,depth)  see ldf_dyn_c3d.h90
55      !!
56      !!      N.B. User defined include files.  By default, 3d and 2d coef.
57      !!      are set to a constant value given in the namelist and the 1d
58      !!      coefficients are initialized to a hyperbolic tangent vertical
59      !!      profile.
60      !!
61      !! Reference :   Madec, G. and M. Imbard, 1996: Climate Dynamics, 12, 381-388.
62      !!----------------------------------------------------------------------
63      INTEGER ::   ioptio         ! ???
64      LOGICAL ::   ll_print = .FALSE.    ! Logical flag for printing viscosity coef.
65      !!
66      NAMELIST/namdyn_ldf/ ln_dynldf_lap  , ln_dynldf_bilap,                  &
67         &                 ln_dynldf_level, ln_dynldf_hor  , ln_dynldf_iso,   &
68         &                 rn_ahm_0_lap   , rn_ahmb_0      , rn_ahm_0_blp
69      !!----------------------------------------------------------------------
70
71      REWIND( numnam )                  ! Read Namelist namdyn_ldf : Lateral physics
72      READ  ( numnam, namdyn_ldf )
73
74      IF(lwp) THEN                      ! Parameter print
75         WRITE(numout,*)
76         WRITE(numout,*) 'ldf_dyn : lateral momentum physics'
77         WRITE(numout,*) '~~~~~~~'
78         WRITE(numout,*) '   Namelist nam_dynldf : set lateral mixing parameters'
79         WRITE(numout,*) '      laplacian operator                      ln_dynldf_lap   = ', ln_dynldf_lap
80         WRITE(numout,*) '      bilaplacian operator                    ln_dynldf_bilap = ', ln_dynldf_bilap
81         WRITE(numout,*) '      iso-level                               ln_dynldf_level = ', ln_dynldf_level
82         WRITE(numout,*) '      horizontal (geopotential)               ln_dynldf_hor   = ', ln_dynldf_hor
83         WRITE(numout,*) '      iso-neutral                             ln_dynldf_iso   = ', ln_dynldf_iso
84         WRITE(numout,*) '      horizontal laplacian eddy viscosity     rn_ahm_0_lap    = ', rn_ahm_0_lap
85         WRITE(numout,*) '      background viscosity                    rn_ahmb_0       = ', rn_ahmb_0
86         WRITE(numout,*) '      horizontal bilaplacian eddy viscosity   rn_ahm_0_blp    = ', rn_ahm_0_blp
87      ENDIF
88
89      ahm0     = rn_ahm_0_lap              ! OLD namelist variables defined from DOCTOR namelist variables
90      ahmb0    = rn_ahmb_0
91      ahm0_blp = rn_ahm_0_blp
92
93      ! ... check of lateral diffusive operator on tracers
94      !           ==> will be done in trazdf module
95
96      ! ... Space variation of eddy coefficients
97      ioptio = 0
98#if defined key_dynldf_c3d
99      IF(lwp) WRITE(numout,*) '   momentum mixing coef. = F( latitude, longitude, depth)'
100      ioptio = ioptio+1
101#endif
102#if defined key_dynldf_c2d
103      IF(lwp) WRITE(numout,*) '   momentum mixing coef. = F( latitude, longitude)'
104      ioptio = ioptio+1
105#endif
106#if defined key_dynldf_c1d
107      IF(lwp) WRITE(numout,*) '   momentum mixing coef. = F( depth )'
108      ioptio = ioptio+1
109      IF( ln_sco ) CALL ctl_stop( 'key_dynldf_c1d cannot be used in s-coordinate (ln_sco)' )
110#endif
111      IF( ioptio == 0 ) THEN
112          IF(lwp) WRITE(numout,*) '   momentum mixing coef. = constant  (default option)'
113        ELSEIF( ioptio > 1 ) THEN
114           CALL ctl_stop( 'use only one of the following keys: key_dynldf_c3d, key_dynldf_c2d, key_dynldf_c1d' )
115      ENDIF
116
117
118      IF( ln_dynldf_bilap ) THEN
119         IF(lwp) WRITE(numout,*) '   biharmonic momentum diffusion'
120         IF( .NOT. ln_dynldf_lap ) ahm0 = ahm0_blp   ! Allow spatially varying coefs, which use ahm0 as input
121         IF( ahm0_blp > 0 .AND. .NOT. lk_esopa )   CALL ctl_stop( 'The horizontal viscosity coef. ahm0 must be negative' )
122      ELSE
123         IF(lwp) WRITE(numout,*) '   harmonic momentum diff. (default)'
124         IF( ahm0 < 0 .AND. .NOT. lk_esopa )   CALL ctl_stop( 'The horizontal viscosity coef. ahm0 must be positive' )
125      ENDIF
126
127
128      ! Lateral eddy viscosity
129      ! ======================
130#if defined key_dynldf_c3d
131      CALL ldf_dyn_c3d( ll_print )   ! ahm = 3D coef. = F( longitude, latitude, depth )
132#elif defined key_dynldf_c2d
133      CALL ldf_dyn_c2d( ll_print )   ! ahm = 1D coef. = F( longitude, latitude )
134#elif defined key_dynldf_c1d
135      CALL ldf_dyn_c1d( ll_print )   ! ahm = 1D coef. = F( depth )
136#else
137                                     ! Constant coefficients
138      IF(lwp) WRITE(numout,*)
139      IF(lwp) WRITE(numout,*) 'inildf: constant eddy viscosity coef. '
140      IF(lwp) WRITE(numout,*) '~~~~~~'
141      IF(lwp) WRITE(numout,*) '        ahm1 = ahm2 = ahm0 =  ',ahm0
142#endif
143      !
144   END SUBROUTINE ldf_dyn_init
145
146#if defined key_dynldf_c3d
147#   include "ldfdyn_c3d.h90"
148#elif defined key_dynldf_c2d
149#   include "ldfdyn_c2d.h90"
150#elif defined key_dynldf_c1d
151#   include "ldfdyn_c1d.h90"
152#endif
153
154
155   SUBROUTINE ldf_zpf_1d( ld_print, pdam, pwam, pbot, pdep, pah )
156      !!----------------------------------------------------------------------
157      !!                  ***  ROUTINE ldf_zpf  ***
158      !!                   
159      !! ** Purpose :   vertical adimensional profile for eddy coefficient
160      !!
161      !! ** Method  :   1D eddy viscosity coefficients ( depth )
162      !!----------------------------------------------------------------------
163      LOGICAL , INTENT(in   )                 ::   ld_print   ! If true, output arrays on numout
164      REAL(wp), INTENT(in   )                 ::   pdam       ! depth of the inflection point
165      REAL(wp), INTENT(in   )                 ::   pwam       ! width of inflection
166      REAL(wp), INTENT(in   )                 ::   pbot       ! bottom value (0<pbot<= 1)
167      REAL(wp), INTENT(in   ), DIMENSION(jpk) ::   pdep       ! depth of the gridpoint (T, U, V, F)
168      REAL(wp), INTENT(inout), DIMENSION(jpk) ::   pah        ! adimensional vertical profile
169      !!
170      INTEGER  ::   jk           ! dummy loop indices
171      REAL(wp) ::   zm00, zm01, zmhb, zmhs       ! temporary scalars
172      !!----------------------------------------------------------------------
173
174      zm00 = TANH( ( pdam - gdept_0(1    ) ) / pwam )
175      zm01 = TANH( ( pdam - gdept_0(jpkm1) ) / pwam )
176      zmhs = zm00 / zm01
177      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
178
179      DO jk = 1, jpk
180         pah(jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam )  )
181      END DO
182
183      IF(lwp .AND. ld_print ) THEN      ! Control print
184         WRITE(numout,*)
185         WRITE(numout,*) '         ahm profile : '
186         WRITE(numout,*)
187         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
188         DO jk = 1, jpk
189            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(jk), pdep(jk)
190         END DO
191      ENDIF
192      !
193   END SUBROUTINE ldf_zpf_1d
194
195
196   SUBROUTINE ldf_zpf_1d_3d( ld_print, pdam, pwam, pbot, pdep, pah )
197      !!----------------------------------------------------------------------
198      !!                  ***  ROUTINE ldf_zpf  ***
199      !!
200      !! ** Purpose :   vertical adimensional profile for eddy coefficient
201      !!
202      !! ** Method  :   1D eddy viscosity coefficients ( depth )
203      !!----------------------------------------------------------------------
204      LOGICAL , INTENT(in   )                         ::   ld_print   ! If true, output arrays on numout
205      REAL(wp), INTENT(in   )                         ::   pdam       ! depth of the inflection point
206      REAL(wp), INTENT(in   )                         ::   pwam       ! width of inflection
207      REAL(wp), INTENT(in   )                         ::   pbot       ! bottom value (0<pbot<= 1)
208      REAL(wp), INTENT(in   ), DIMENSION          (:) ::   pdep       ! depth of the gridpoint (T, U, V, F)
209      REAL(wp), INTENT(inout), DIMENSION      (:,:,:) ::   pah        ! adimensional vertical profile
210      !!
211      INTEGER  ::   jk           ! dummy loop indices
212      REAL(wp) ::   zm00, zm01, zmhb, zmhs, zcf  ! temporary scalars
213      !!----------------------------------------------------------------------
214
215      zm00 = TANH( ( pdam - gdept_0(1    ) ) / pwam )
216      zm01 = TANH( ( pdam - gdept_0(jpkm1) ) / pwam )
217      zmhs = zm00 / zm01
218      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
219
220      DO jk = 1, jpk
221         zcf = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam )  )
222         pah(:,:,jk) = zcf
223      END DO
224
225      IF(lwp .AND. ld_print ) THEN      ! Control print
226         WRITE(numout,*)
227         WRITE(numout,*) '         ahm profile : '
228         WRITE(numout,*)
229         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
230         DO jk = 1, jpk
231            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(jk)
232         END DO
233      ENDIF
234      !
235   END SUBROUTINE ldf_zpf_1d_3d
236
237
238   SUBROUTINE ldf_zpf_3d( ld_print, pdam, pwam, pbot, pdep, pah )
239      !!----------------------------------------------------------------------
240      !!                  ***  ROUTINE ldf_zpf  ***
241      !!
242      !! ** Purpose :   vertical adimensional profile for eddy coefficient
243      !!
244      !! ** Method  :   3D for partial step or s-coordinate
245      !!----------------------------------------------------------------------
246      LOGICAL , INTENT(in   )                         ::   ld_print   ! If true, output arrays on numout
247      REAL(wp), INTENT(in   )                         ::   pdam       ! depth of the inflection point
248      REAL(wp), INTENT(in   )                         ::   pwam       ! width of inflection
249      REAL(wp), INTENT(in   )                         ::   pbot       ! bottom value (0<pbot<= 1)
250      REAL(wp), INTENT(in   ), DIMENSION      (:,:,:) ::   pdep       ! dep of the gridpoint (T, U, V, F)
251      REAL(wp), INTENT(inout), DIMENSION      (:,:,:) ::   pah        ! adimensional vertical profile
252      !!
253      INTEGER  ::   jk           ! dummy loop indices
254      REAL(wp) ::   zm00, zm01, zmhb, zmhs       ! temporary scalars
255      !!----------------------------------------------------------------------
256
257      zm00 = TANH( ( pdam - gdept_0(1    ) ) / pwam )   
258      zm01 = TANH( ( pdam - gdept_0(jpkm1) ) / pwam )
259      zmhs = zm00 / zm01
260      zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01
261
262      DO jk = 1, jpk
263         pah(:,:,jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(:,:,jk) ) / pwam )  )
264      END DO
265
266      IF(lwp .AND. ld_print ) THEN      ! Control print
267         WRITE(numout,*)
268         WRITE(numout,*) '         ahm profile : '
269         WRITE(numout,*)
270         WRITE(numout,'("  jk      ahm       ","  depth t-level " )')
271         DO jk = 1, jpk
272            WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(1,1,jk)
273         END DO
274      ENDIF
275      !
276   END SUBROUTINE ldf_zpf_3d
277
278   !!======================================================================
279END MODULE ldfdyn
Note: See TracBrowser for help on using the repository browser.