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

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

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • 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(jpk) ::   pdep       ! depth of the gridpoint (T, U, V, F)
174      REAL(wp), INTENT(inout), DIMENSION(jpk) ::   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.