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/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/LDF/ldfc1d_c2d.F90 @ 7037

Last change on this file since 7037 was 7037, checked in by mocavero, 8 years ago

ORCA2_LIM_PISCES hybrid version update
File size: 10.9 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!$OMP PARALLEL DO schedule(static) private(jj,ji,zemax1,zemax2)
144            DO jj = 1, jpj 
145               DO ji = 1, jpi 
146                  zemax1 = MAX( e1t(ji,jj), e2t(ji,jj) ) * tmask(ji,jj,1)
147                  zemax2 = MAX( e1f(ji,jj), e2f(ji,jj) ) * fmask(ji,jj,1)
148                  pah1(ji,jj,1) = za00 * zemax1
149                  pah2(ji,jj,1) = za00 * zemax2
150               END DO
151            END DO
152         ELSEIF( cd_op == 'BLP' ) THEN     ! bilaplacian operator
153            IF(lwp) WRITE(numout,*) '              momentum bilaplacian coeffcients = rn_bht0/e_equ * max(e1,e2)**3'
154            za00 = pah0 / ( zd_max * zd_max * zd_max )
155!$OMP PARALLEL DO schedule(static) private(jj,ji,zemax1,zemax2)
156            DO jj = 1, jpj
157               DO ji = 1, jpi
158                  zemax1 = MAX( e1t(ji,jj), e2t(ji,jj) ) * tmask(ji,jj,1)
159                  zemax2 = MAX( e1f(ji,jj), e2f(ji,jj) ) * fmask(ji,jj,1)
160                  pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1 
161                  pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2 
162               END DO
163            END DO
164         ELSE                                ! NO diffusion/viscosity
165            CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' )
166         ENDIF
167         !                                !  deeper values  (LAP and BLP cases)
168!$OMP PARALLEL DO schedule(static) private(jk)
169         DO jk = 2, jpk
170            pah1(:,:,jk) = pah1(:,:,1) * tmask(:,:,jk) 
171            pah2(:,:,jk) = pah2(:,:,1) * fmask(:,:,jk) 
172         END DO
173         !
174      CASE( 'TRA' )                     ! U- and V-points (approximation in zps-coord.)
175         IF( cd_op == 'LAP' ) THEN            ! laplacian operator
176            IF(lwp) WRITE(numout,*) '              tracer laplacian coeffcients = rn_aht0/e_equ * max(e1,e2)'
177            za00 = pah0 / zd_max
178!$OMP PARALLEL DO schedule(static) private(jj,ji,zemax1,zemax2)
179            DO jj = 1, jpj 
180               DO ji = 1, jpi 
181                  zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1)
182                  zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1)
183                  pah1(ji,jj,1) = za00 * zemax1
184                  pah2(ji,jj,1) = za00 * zemax2
185               END DO
186            END DO
187         ELSEIF( cd_op == 'BLP' ) THEN      ! bilaplacian operator (NB: square root of the coeff)
188            IF(lwp) WRITE(numout,*) '              tracer bilaplacian coeffcients = rn_bht0/e_equ * max(e1,e2)**3'
189            za00 = pah0 / ( zd_max * zd_max * zd_max )
190!$OMP PARALLEL DO schedule(static) private(jj,ji,zemax1,zemax2)
191            DO jj = 1, jpj
192               DO ji = 1, jpi
193                  zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1) 
194                  zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1) 
195                  pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1 
196                  pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2 
197               END DO
198            END DO
199         ELSE                                ! NO diffusion/viscosity
200            CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' )
201         ENDIF
202         !                                !  deeper values  (LAP and BLP cases)
203!$OMP PARALLEL DO schedule(static) private(jk)
204         DO jk = 2, jpk
205            pah1(:,:,jk) = pah1(:,:,1) * umask(:,:,jk) 
206            pah2(:,:,jk) = pah2(:,:,1) * vmask(:,:,jk) 
207         END DO
208         !
209      END SELECT
210      !
211   END SUBROUTINE ldf_c2d
212
213   !!======================================================================
214END MODULE ldfc1d_c2d
Note: See TracBrowser for help on using the repository browser.