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/2017/dev_r8126_ROBUST08_no_ghost/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/2017/dev_r8126_ROBUST08_no_ghost/NEMOGCM/NEMO/OPA_SRC/LDF/ldfc1d_c2d.F90 @ 8809

Last change on this file since 8809 was 7753, checked in by mocavero, 7 years ago

Reverting trunk to remove OpenMP

File size: 10.8 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      IF(lwp) THEN
63         WRITE(numout,*)
64         WRITE(numout,*) '   ldf_c1d : set a given profile to eddy diffusivity/viscosity coefficients'
65         WRITE(numout,*) '   ~~~~~~~'
66      ENDIF
67
68      ! initialization of the profile
69      zh =  500._wp              ! depth    of the inflection point [m]
70      zw =  1._wp / 200._wp      ! width^-1     -        -      -   [1/m]
71      !                          ! associated coefficient           [-]
72      zc = ( 1._wp - prat ) / ( 1._wp + TANH( zh * zw) )
73      !
74      !
75      SELECT CASE( cd_type )        ! point of calculation
76      !
77      CASE( 'DYN' )                     ! T- and F-points
78         DO jk = 1, jpk                      ! pah1 at T-point
79            pah1(:,:,jk) = pahs1(:,:) * (  prat + zc * ( 1._wp + TANH( - ( gdept_n(:,:,jk) - zh ) * zw) )  ) * tmask(:,:,jk)
80         END DO
81         DO jk = 1, jpk                      ! pah2 at F-point (zdep2 is an approximation in zps-coord.)
82            DO jj = 1, jpjm1
83               DO ji = 1, fs_jpim1
84                  zdep2 = (  gdept_n(ji,jj+1,jk) + gdept_n(ji+1,jj+1,jk)   &
85                     &     + gdept_n(ji,jj  ,jk) + gdept_n(ji+1,jj  ,jk)  ) * 0.25_wp
86                  pah2(ji,jj,jk) = pahs2(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  ) * fmask(ji,jj,jk)
87               END DO
88            END DO
89         END DO
90         CALL lbc_lnk( pah2, 'F', 1. )   ! Lateral boundary conditions
91         !
92      CASE( 'TRA' )                     ! U- and V-points (zdep1 & 2 are an approximation in zps-coord.)
93         DO jk = 1, jpk
94            DO jj = 1, jpjm1
95               DO ji = 1, fs_jpim1
96                  zdep1 = (  gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk)  ) * 0.5_wp
97                  zdep2 = (  gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk)  ) * 0.5_wp
98                  pah1(ji,jj,jk) = pahs1(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) )  ) * umask(ji,jj,jk)
99                  pah2(ji,jj,jk) = pahs2(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  ) * vmask(ji,jj,jk)
100               END DO
101            END DO
102         END DO
103         CALL lbc_lnk( pah1, 'U', 1. )   ! Lateral boundary conditions
104         CALL lbc_lnk( pah2, 'V', 1. )   
105         !
106      END SELECT
107      !
108   END SUBROUTINE ldf_c1d
109
110
111   SUBROUTINE ldf_c2d( cd_type, cd_op, pah0, pah1, pah2 )
112      !!----------------------------------------------------------------------
113      !!                  ***  ROUTINE ldf_c2d  ***
114      !!             
115      !! ** Purpose :   2D eddy diffusivity/viscosity coefficients
116      !!
117      !! ** Method  :   2D eddy diffusivity coefficients F( e1 , e2 )
118      !!       laplacian   operator :   ah proportional to the scale factor      [m2/s]
119      !!       bilaplacian operator :   ah proportional to the (scale factor)^3  [m4/s]
120      !!       In both cases, pah0 is the maximum value reached by the coefficient
121      !!       at the Equator in case of e1=ra*rad= ~111km, not over the whole domain.
122      !!
123      !!   cd_type = TRA      pah1, pah2 defined at U- and V-points
124      !!             DYN      pah1, pah2 defined at T- and F-points
125      !!----------------------------------------------------------------------
126      CHARACTER(len=3)                , INTENT(in   ) ::   cd_type      ! DYNamique or TRAcers
127      CHARACTER(len=3)                , INTENT(in   ) ::   cd_op        ! operator: LAPlacian BiLaPlacian
128      REAL(wp)                        , INTENT(in   ) ::   pah0         ! eddy coefficient   [m2/s or m4/s]
129      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pah1, pah2   ! eddy coefficient   [m2/s or m4/s]
130      !
131      INTEGER  ::   ji, jj, jk   ! dummy loop indices
132      REAL(wp) ::   za00, zd_max, zemax1, zemax2   ! local scalar
133      !!----------------------------------------------------------------------
134      !
135      zd_max = ra * rad       ! = 1 degree at the equator in meters
136      !
137      IF(lwp) THEN
138         WRITE(numout,*)
139         WRITE(numout,*) '   ldf_c2d :     aht = rn_aht0 *  max(e1,e2)/e_equator     (  laplacian) '
140         WRITE(numout,*) '   ~~~~~~~       or  = rn_bht0 * [max(e1,e2)/e_equator]**3 (bilaplacian)'
141         WRITE(numout,*)
142      ENDIF
143      !
144      SELECT CASE( cd_type )        !==  surface values  ==!  (depending on DYN/TRA)
145      !
146      CASE( 'DYN' )                     ! T- and F-points
147         IF( cd_op == 'LAP' ) THEN            ! laplacian operator
148            IF(lwp) WRITE(numout,*) '              momentum laplacian coeffcients = rn_aht0/e_equ * max(e1,e2)'
149            za00 = pah0 / zd_max
150            DO jj = 1, jpj 
151               DO ji = 1, jpi 
152                  zemax1 = MAX( e1t(ji,jj), e2t(ji,jj) ) * tmask(ji,jj,1)
153                  zemax2 = MAX( e1f(ji,jj), e2f(ji,jj) ) * fmask(ji,jj,1)
154                  pah1(ji,jj,1) = za00 * zemax1
155                  pah2(ji,jj,1) = za00 * zemax2
156               END DO
157            END DO
158         ELSEIF( cd_op == 'BLP' ) THEN     ! bilaplacian operator
159            IF(lwp) WRITE(numout,*) '              momentum bilaplacian coeffcients = rn_bht0/e_equ * max(e1,e2)**3'
160            za00 = pah0 / ( zd_max * zd_max * zd_max )
161            DO jj = 1, jpj
162               DO ji = 1, jpi
163                  zemax1 = MAX( e1t(ji,jj), e2t(ji,jj) ) * tmask(ji,jj,1)
164                  zemax2 = MAX( e1f(ji,jj), e2f(ji,jj) ) * fmask(ji,jj,1)
165                  pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1 
166                  pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2 
167               END DO
168            END DO
169         ELSE                                ! NO diffusion/viscosity
170            CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' )
171         ENDIF
172         !                                !  deeper values  (LAP and BLP cases)
173         DO jk = 2, jpk
174            pah1(:,:,jk) = pah1(:,:,1) * tmask(:,:,jk) 
175            pah2(:,:,jk) = pah2(:,:,1) * fmask(:,:,jk) 
176         END DO
177         !
178      CASE( 'TRA' )                     ! U- and V-points (approximation in zps-coord.)
179         IF( cd_op == 'LAP' ) THEN            ! laplacian operator
180            IF(lwp) WRITE(numout,*) '              tracer laplacian coeffcients = rn_aht0/e_equ * max(e1,e2)'
181            za00 = pah0 / zd_max
182            DO jj = 1, jpj 
183               DO ji = 1, jpi 
184                  zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1)
185                  zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1)
186                  pah1(ji,jj,1) = za00 * zemax1
187                  pah2(ji,jj,1) = za00 * zemax2
188               END DO
189            END DO
190         ELSEIF( cd_op == 'BLP' ) THEN      ! bilaplacian operator (NB: square root of the coeff)
191            IF(lwp) WRITE(numout,*) '              tracer bilaplacian coeffcients = rn_bht0/e_equ * max(e1,e2)**3'
192            za00 = pah0 / ( zd_max * zd_max * zd_max )
193            DO jj = 1, jpj
194               DO ji = 1, jpi
195                  zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1) 
196                  zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1) 
197                  pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1 
198                  pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2 
199               END DO
200            END DO
201         ELSE                                ! NO diffusion/viscosity
202            CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' )
203         ENDIF
204         !                                !  deeper values  (LAP and BLP cases)
205         DO jk = 2, jpk
206            pah1(:,:,jk) = pah1(:,:,1) * umask(:,:,jk) 
207            pah2(:,:,jk) = pah2(:,:,1) * vmask(:,:,jk) 
208         END DO
209         !
210      END SELECT
211      !
212   END SUBROUTINE ldf_c2d
213
214   !!======================================================================
215END MODULE ldfc1d_c2d
Note: See TracBrowser for help on using the repository browser.