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

source: branches/2012/dev_r3452_NOCL02_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 3497

Last change on this file since 3497 was 3497, checked in by hliu, 12 years ago

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