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

source: branches/DEV_r1784_mid_year_merge_2010/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 1954

Last change on this file since 1954 was 1954, checked in by acc, 14 years ago

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