source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 4409

Last change on this file since 4409 was 4409, checked in by trackstand2, 7 years ago

Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.

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