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.
ldfc1d_c2d.F90 in branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/LDF/ldfc1d_c2d.F90 @ 5845

Last change on this file since 5845 was 5845, checked in by gm, 8 years ago

#1613: vvl by default: suppression of domzgr_substitute.h90

File size: 10.6 KB
Line 
1MODULE ldfc1d_c2d
2   !!======================================================================
3   !!                    ***  MODULE  ldfc1d_c2d  ***
4   !! Ocean physics:  profile and horizontal shape of lateral eddy coefficients
5   !!=====================================================================
6   !! History :  3.7  ! 2013-12  (G. Madec)  restructuration/simplification of aht/aeiv specification,
7   !!                 !                      add velocity dependent coefficient and optional read in file
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   ldf_c1d       : ah reduced by 1/4 on the vertical (tanh profile, inflection at 300m)
12   !!   ldf_c2d       : ah = F(e1,e2) (laplacian or = F(e1^3,e2^3) (bilaplacian)
13   !!----------------------------------------------------------------------
14   USE oce            ! ocean dynamics and tracers
15   USE dom_oce        ! ocean space and time domain
16   USE phycst         ! physical constants
17   !
18   USE in_out_manager ! I/O manager
19   USE lib_mpp        ! distribued memory computing library
20   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
21
22   IMPLICIT NONE
23   PRIVATE
24
25   PUBLIC   ldf_c1d   ! called by ldftra and ldfdyn modules
26   PUBLIC   ldf_c2d   ! called by ldftra and ldfdyn modules
27
28 
29   !! * Substitutions
30#  include "vectopt_loop_substitute.h90"
31   !!----------------------------------------------------------------------
32   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
33   !! $Id: $
34   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
35   !!----------------------------------------------------------------------
36CONTAINS
37
38   SUBROUTINE ldf_c1d( cd_type, prat, pahs1, pahs2, pah1, pah2 )
39      !!----------------------------------------------------------------------
40      !!                  ***  ROUTINE ldf_c1d  ***
41      !!             
42      !! ** Purpose :   1D eddy diffusivity/viscosity coefficients
43      !!
44      !! ** Method  :   1D eddy diffusivity coefficients F( depth )
45      !!                Reduction by prat from surface to bottom
46      !!                hyperbolic tangent profile with inflection point
47      !!                at zh=500m and a width of zw=200m
48      !!
49      !!   cd_type = TRA      pah1, pah2 defined at U- and V-points
50      !!             DYN      pah1, pah2 defined at T- and F-points
51      !!----------------------------------------------------------------------
52      CHARACTER(len=2)                , INTENT(in   ) ::   cd_type        ! DYNamique or TRAcers
53      REAL(wp)                        , INTENT(in   ) ::   prat           ! ratio surface/deep values           [-]
54      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   pahs1, pahs2   ! surface value of eddy coefficient   [m2/s]
55      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pah1 , pah2    ! eddy coefficient                    [m2/s]
56      !
57      INTEGER  ::   ji, jj, jk      ! dummy loop indices
58      REAL(wp) ::   zh, zc, zdep1   ! local scalars
59      REAL(wp) ::   zw    , zdep2   !   -      -
60      !!----------------------------------------------------------------------
61
62      ! initialization of the profile
63      zh =  500._wp              ! depth    of the inflection point [m]
64      zw =  1._wp / 200._wp      ! width^-1     -        -      -   [1/m]
65      !                          ! associated coefficient           [-]
66      zc = ( 1._wp - prat ) / ( 1._wp + TANH( zh * zw) )
67      !
68      !
69      SELECT CASE( cd_type )        ! point of calculation
70      !
71      CASE( 'DYN' )                     ! T- and F-points
72         DO jk = 1, jpk                      ! pah1 at T-point
73            pah1(:,:,jk) = pahs1(:,:) * (  prat + zc * ( 1._wp + TANH( - ( gdept_n(:,:,jk) - zh ) * zw) )  ) * tmask(:,:,jk)
74         END DO
75         DO jk = 1, jpk                      ! pah2 at F-point (zdep2 is an approximation in zps-coord.)
76            DO jj = 1, jpjm1
77               DO ji = 1, fs_jpim1
78                  zdep2 = (  gdept_n(ji,jj+1,jk) + gdept_n(ji+1,jj+1,jk)   &
79                     &     + gdept_n(ji,jj  ,jk) + gdept_n(ji+1,jj  ,jk)  ) * 0.25_wp
80                  pah2(ji,jj,jk) = pahs2(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  ) * fmask(ji,jj,jk)
81               END DO
82            END DO
83         END DO
84         CALL lbc_lnk( pah2, 'F', 1. )   ! Lateral boundary conditions
85         !
86      CASE( 'TRA' )                     ! U- and V-points (zdep1 & 2 are an approximation in zps-coord.)
87         DO jk = 1, jpk
88            DO jj = 1, jpjm1
89               DO ji = 1, fs_jpim1
90                  zdep1 = (  gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk)  ) * 0.5_wp
91                  zdep2 = (  gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk)  ) * 0.5_wp
92                  pah1(ji,jj,jk) = pahs1(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) )  ) * umask(ji,jj,jk)
93                  pah2(ji,jj,jk) = pahs2(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  ) * vmask(ji,jj,jk)
94               END DO
95            END DO
96         END DO
97         CALL lbc_lnk( pah1, 'U', 1. )   ! Lateral boundary conditions
98         CALL lbc_lnk( pah2, 'V', 1. )   
99         !
100      END SELECT
101      !
102   END SUBROUTINE ldf_c1d
103
104
105   SUBROUTINE ldf_c2d( cd_type, cd_op, pah0, pah1, pah2 )
106      !!----------------------------------------------------------------------
107      !!                  ***  ROUTINE ldf_c2d  ***
108      !!             
109      !! ** Purpose :   2D eddy diffusivity/viscosity coefficients
110      !!
111      !! ** Method  :   2D eddy diffusivity coefficients F( e1 , e2 )
112      !!       laplacian   operator :   ah proportional to the scale factor      [m2/s]
113      !!       bilaplacian operator :   ah proportional to the (scale factor)^3  [m4/s]
114      !!       In both cases, pah0 is the maximum value reached by the coefficient
115      !!       at the Equator in case of e1=ra*rad= ~111km, not over the whole domain.
116      !!
117      !!   cd_type = TRA      pah1, pah2 defined at U- and V-points
118      !!             DYN      pah1, pah2 defined at T- and F-points
119      !!----------------------------------------------------------------------
120      CHARACTER(len=3)                , INTENT(in   ) ::   cd_type      ! DYNamique or TRAcers
121      CHARACTER(len=3)                , INTENT(in   ) ::   cd_op        ! operator: LAPlacian BiLaPlacian
122      REAL(wp)                        , INTENT(in   ) ::   pah0         ! eddy coefficient   [m2/s or m4/s]
123      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pah1, pah2   ! eddy coefficient   [m2/s or m4/s]
124      !
125      INTEGER  ::   ji, jj, jk   ! dummy loop indices
126      REAL(wp) ::   za00, zd_max, zemax1, zemax2   ! local scalar
127      !!----------------------------------------------------------------------
128      !
129      zd_max = ra * rad       ! = 1 degree at the equator in meters
130      !
131      IF(lwp) THEN
132         WRITE(numout,*) 'ldf_c2d :     aht = rn_aht0 *  max(e1,e2)/e_equator     (  laplacian) '
133         WRITE(numout,*) '~~~~~~~       or  = rn_bht0 * [max(e1,e2)/e_equator]**3 (bilaplacian)'
134         WRITE(numout,*)
135      ENDIF
136      !
137      SELECT CASE( cd_type )        !==  surface values  ==!  (depending on DYN/TRA)
138      !
139      CASE( 'DYN' )                     ! T- and F-points
140         IF( cd_op == 'LAP' ) THEN            ! laplacian operator
141            IF(lwp) WRITE(numout,*) '              momentum laplacian coeffcients = rn_aht0/e_equ * max(e1,e2)'
142            za00 = pah0 / zd_max
143            DO jj = 1, jpj 
144               DO ji = 1, jpi 
145                  zemax1 = MAX( e1t(ji,jj), e2t(ji,jj) ) * tmask(ji,jj,1)
146                  zemax2 = MAX( e1f(ji,jj), e2f(ji,jj) ) * fmask(ji,jj,1)
147                  pah1(ji,jj,1) = za00 * zemax1
148                  pah2(ji,jj,1) = za00 * zemax2
149               END DO
150            END DO
151         ELSEIF( cd_op == 'BLP' ) THEN     ! bilaplacian operator
152            IF(lwp) WRITE(numout,*) '              momentum bilaplacian coeffcients = rn_bht0/e_equ * max(e1,e2)**3'
153            za00 = pah0 / ( zd_max * zd_max * zd_max )
154            DO jj = 1, jpj
155               DO ji = 1, jpi
156                  zemax1 = MAX( e1t(ji,jj), e2t(ji,jj) ) * tmask(ji,jj,1)
157                  zemax2 = MAX( e1f(ji,jj), e2f(ji,jj) ) * fmask(ji,jj,1)
158                  pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1 
159                  pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2 
160               END DO
161            END DO
162         ELSE                                ! NO diffusion/viscosity
163            CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' )
164         ENDIF
165         !                                !  deeper values  (LAP and BLP cases)
166         DO jk = 2, jpk
167            pah1(:,:,jk) = pah1(:,:,1) * tmask(:,:,jk) 
168            pah2(:,:,jk) = pah2(:,:,1) * fmask(:,:,jk) 
169         END DO
170         !
171      CASE( 'TRA' )                     ! U- and V-points (approximation in zps-coord.)
172         IF( cd_op == 'LAP' ) THEN            ! laplacian operator
173            IF(lwp) WRITE(numout,*) '              tracer laplacian coeffcients = rn_aht0/e_equ * max(e1,e2)'
174            za00 = pah0 / zd_max
175            DO jj = 1, jpj 
176               DO ji = 1, jpi 
177                  zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1)
178                  zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1)
179                  pah1(ji,jj,1) = za00 * zemax1
180                  pah2(ji,jj,1) = za00 * zemax2
181               END DO
182            END DO
183         ELSEIF( cd_op == 'BLP' ) THEN      ! bilaplacian operator (NB: square root of the coeff)
184            IF(lwp) WRITE(numout,*) '              tracer bilaplacian coeffcients = rn_bht0/e_equ * max(e1,e2)**3'
185            za00 = pah0 / ( zd_max * zd_max * zd_max )
186            DO jj = 1, jpj
187               DO ji = 1, jpi
188                  zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1) 
189                  zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1) 
190                  pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1 
191                  pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2 
192               END DO
193            END DO
194         ELSE                                ! NO diffusion/viscosity
195            CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' )
196         ENDIF
197         !                                !  deeper values  (LAP and BLP cases)
198         DO jk = 2, jpk
199            pah1(:,:,jk) = pah1(:,:,1) * umask(:,:,jk) 
200            pah2(:,:,jk) = pah2(:,:,1) * vmask(:,:,jk) 
201         END DO
202         !
203      END SELECT
204      !
205   END SUBROUTINE ldf_c2d
206
207   !!======================================================================
208END MODULE ldfc1d_c2d
Note: See TracBrowser for help on using the repository browser.