source: trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfc1d_c2d.F90 @ 5836

Last change on this file since 5836 was 5836, checked in by cetlod, 5 years ago

merge the simplification branch onto the trunk, see ticket #1612

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 "domzgr_substitute.h90"
31#  include "vectopt_loop_substitute.h90"
32   !!----------------------------------------------------------------------
33   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
34   !! $Id: $
35   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
36   !!----------------------------------------------------------------------
37CONTAINS
38
39   SUBROUTINE ldf_c1d( cd_type, prat, pahs1, pahs2, pah1, pah2 )
40      !!----------------------------------------------------------------------
41      !!                  ***  ROUTINE ldf_c1d  ***
42      !!             
43      !! ** Purpose :   1D eddy diffusivity/viscosity coefficients
44      !!
45      !! ** Method  :   1D eddy diffusivity coefficients F( depth )
46      !!                Reduction by prat from surface to bottom
47      !!                hyperbolic tangent profile with inflection point
48      !!                at zh=500m and a width of zw=200m
49      !!
50      !!   cd_type = TRA      pah1, pah2 defined at U- and V-points
51      !!             DYN      pah1, pah2 defined at T- and F-points
52      !!----------------------------------------------------------------------
53      CHARACTER(len=2)                , INTENT(in   ) ::   cd_type        ! DYNamique or TRAcers
54      REAL(wp)                        , INTENT(in   ) ::   prat           ! ratio surface/deep values           [-]
55      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   pahs1, pahs2   ! surface value of eddy coefficient   [m2/s]
56      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pah1 , pah2    ! eddy coefficient                    [m2/s]
57      !
58      INTEGER  ::   ji, jj, jk      ! dummy loop indices
59      REAL(wp) ::   zh, zc, zdep1   ! local scalars
60      REAL(wp) ::   zw    , zdep2   !   -      -
61      !!----------------------------------------------------------------------
62
63      ! initialization of the profile
64      zh =  500._wp              ! depth    of the inflection point [m]
65      zw =  1._wp / 200._wp      ! width^-1     -        -      -   [1/m]
66      !                          ! associated coefficient           [-]
67      zc = ( 1._wp - prat ) / ( 1._wp + TANH( zh * zw) )
68      !
69      !
70      SELECT CASE( cd_type )        ! point of calculation
71      !
72      CASE( 'DYN' )                     ! T- and F-points
73         DO jk = 1, jpk                      ! pah1 at T-point
74            pah1(:,:,jk) = pahs1(:,:) * (  prat + zc * ( 1._wp + TANH( - ( fsdept(:,:,jk) - zh ) * zw) )  ) * tmask(:,:,jk)
75         END DO
76         DO jk = 1, jpk                      ! pah2 at F-point (zdep2 is an approximation in zps-coord.)
77            DO jj = 1, jpjm1
78               DO ji = 1, fs_jpim1
79                  zdep2 = (  fsdept(ji,jj+1,jk) + fsdept(ji+1,jj+1,jk)   &
80                     &     + fsdept(ji,jj  ,jk) + fsdept(ji+1,jj  ,jk)  ) * 0.25_wp
81                  pah2(ji,jj,jk) = pahs2(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  ) * fmask(ji,jj,jk)
82               END DO
83            END DO
84         END DO
85         CALL lbc_lnk( pah2, 'F', 1. )   ! Lateral boundary conditions
86         !
87      CASE( 'TRA' )                     ! U- and V-points (zdep1 & 2 are an approximation in zps-coord.)
88         DO jk = 1, jpk
89            DO jj = 1, jpjm1
90               DO ji = 1, fs_jpim1
91                  zdep1 = (  fsdept(ji,jj,jk) + fsdept(ji+1,jj,jk)  ) * 0.5_wp
92                  zdep2 = (  fsdept(ji,jj,jk) + fsdept(ji,jj+1,jk)  ) * 0.5_wp
93                  pah1(ji,jj,jk) = pahs1(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) )  ) * umask(ji,jj,jk)
94                  pah2(ji,jj,jk) = pahs2(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  ) * vmask(ji,jj,jk)
95               END DO
96            END DO
97         END DO
98         CALL lbc_lnk( pah1, 'U', 1. )   ! Lateral boundary conditions
99         CALL lbc_lnk( pah2, 'V', 1. )   
100         !
101      END SELECT
102      !
103   END SUBROUTINE ldf_c1d
104
105
106   SUBROUTINE ldf_c2d( cd_type, cd_op, pah0, pah1, pah2 )
107      !!----------------------------------------------------------------------
108      !!                  ***  ROUTINE ldf_c2d  ***
109      !!             
110      !! ** Purpose :   2D eddy diffusivity/viscosity coefficients
111      !!
112      !! ** Method  :   2D eddy diffusivity coefficients F( e1 , e2 )
113      !!       laplacian   operator :   ah proportional to the scale factor      [m2/s]
114      !!       bilaplacian operator :   ah proportional to the (scale factor)^3  [m4/s]
115      !!       In both cases, pah0 is the maximum value reached by the coefficient
116      !!       at the Equator in case of e1=ra*rad= ~111km, not over the whole domain.
117      !!
118      !!   cd_type = TRA      pah1, pah2 defined at U- and V-points
119      !!             DYN      pah1, pah2 defined at T- and F-points
120      !!----------------------------------------------------------------------
121      CHARACTER(len=3)                , INTENT(in   ) ::   cd_type      ! DYNamique or TRAcers
122      CHARACTER(len=3)                , INTENT(in   ) ::   cd_op        ! operator: LAPlacian BiLaPlacian
123      REAL(wp)                        , INTENT(in   ) ::   pah0         ! eddy coefficient   [m2/s or m4/s]
124      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pah1, pah2   ! eddy coefficient   [m2/s or m4/s]
125      !
126      INTEGER  ::   ji, jj, jk   ! dummy loop indices
127      REAL(wp) ::   za00, zd_max, zemax1, zemax2   ! local scalar
128      !!----------------------------------------------------------------------
129      !
130      zd_max = ra * rad       ! = 1 degree at the equator in meters
131      !
132      IF(lwp) THEN
133         WRITE(numout,*) 'ldf_c2d :     aht = rn_aht0 *  max(e1,e2)/e_equator     (  laplacian) '
134         WRITE(numout,*) '~~~~~~~       or  = rn_bht0 * [max(e1,e2)/e_equator]**3 (bilaplacian)'
135         WRITE(numout,*)
136      ENDIF
137      !
138      SELECT CASE( cd_type )        !==  surface values  ==!  (depending on DYN/TRA)
139      !
140      CASE( 'DYN' )                     ! T- and F-points
141         IF( cd_op == 'LAP' ) THEN            ! laplacian operator
142            IF(lwp) WRITE(numout,*) '              momentum laplacian coeffcients = rn_aht0/e_equ * max(e1,e2)'
143            za00 = pah0 / zd_max
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            DO jj = 1, jpj
156               DO ji = 1, jpi
157                  zemax1 = MAX( e1t(ji,jj), e2t(ji,jj) ) * tmask(ji,jj,1)
158                  zemax2 = MAX( e1f(ji,jj), e2f(ji,jj) ) * fmask(ji,jj,1)
159                  pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1 
160                  pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2 
161               END DO
162            END DO
163         ELSE                                ! NO diffusion/viscosity
164            CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' )
165         ENDIF
166         !                                !  deeper values  (LAP and BLP cases)
167         DO jk = 2, jpk
168            pah1(:,:,jk) = pah1(:,:,1) * tmask(:,:,jk) 
169            pah2(:,:,jk) = pah2(:,:,1) * fmask(:,:,jk) 
170         END DO
171         !
172      CASE( 'TRA' )                     ! U- and V-points (approximation in zps-coord.)
173         IF( cd_op == 'LAP' ) THEN            ! laplacian operator
174            IF(lwp) WRITE(numout,*) '              tracer laplacian coeffcients = rn_aht0/e_equ * max(e1,e2)'
175            za00 = pah0 / zd_max
176            DO jj = 1, jpj 
177               DO ji = 1, jpi 
178                  zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1)
179                  zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1)
180                  pah1(ji,jj,1) = za00 * zemax1
181                  pah2(ji,jj,1) = za00 * zemax2
182               END DO
183            END DO
184         ELSEIF( cd_op == 'BLP' ) THEN      ! bilaplacian operator (NB: square root of the coeff)
185            IF(lwp) WRITE(numout,*) '              tracer bilaplacian coeffcients = rn_bht0/e_equ * max(e1,e2)**3'
186            za00 = pah0 / ( zd_max * zd_max * zd_max )
187            DO jj = 1, jpj
188               DO ji = 1, jpi
189                  zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1) 
190                  zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1) 
191                  pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1 
192                  pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2 
193               END DO
194            END DO
195         ELSE                                ! NO diffusion/viscosity
196            CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' )
197         ENDIF
198         !                                !  deeper values  (LAP and BLP cases)
199         DO jk = 2, jpk
200            pah1(:,:,jk) = pah1(:,:,1) * umask(:,:,jk) 
201            pah2(:,:,jk) = pah2(:,:,1) * vmask(:,:,jk) 
202         END DO
203         !
204      END SELECT
205      !
206   END SUBROUTINE ldf_c2d
207
208   !!======================================================================
209END MODULE ldfc1d_c2d
Note: See TracBrowser for help on using the repository browser.