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

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/LDF/ldfc1d_c2d.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
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.