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 NEMO/branches/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/LDF – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/LDF/ldfc1d_c2d.F90 @ 10986

Last change on this file since 10986 was 10986, checked in by andmirek, 5 years ago

GMED 462 add flush
File size: 8.4 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   REAL(wp) ::   r1_2  = 0.5_wp           ! =1/2
29   REAL(wp) ::   r1_4  = 0.25_wp          ! =1/4
30   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! =1/12
31 
32   !!----------------------------------------------------------------------
33   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
34   !! $Id$
35   !! Software governed by the CeCILL license (see ./LICENSE)
36   !!----------------------------------------------------------------------
37CONTAINS
38
39   SUBROUTINE ldf_c1d( cd_type, 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 zratio 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=3)                , INTENT(in   ) ::   cd_type        ! DYNamique or TRAcers
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      REAL(wp) ::   zratio          !   -      -
61      !!----------------------------------------------------------------------
62      !
63      IF(lwp) THEN
64         WRITE(numout,*)
65         WRITE(numout,*) '   ldf_c1d : set a given profile to eddy mixing coefficients'
66         IF(lflush) CALL FLUSH(numout)
67      ENDIF
68      !
69      ! initialization of the profile
70      zratio = 0.25_wp           ! surface/bottom ratio
71      zh =  500._wp              ! depth    of the inflection point [m]
72      zw =  1._wp / 200._wp      ! width^-1     -        -      -   [1/m]
73      !                          ! associated coefficient           [-]
74      zc = ( 1._wp - zratio ) / ( 1._wp + TANH( zh * zw) )
75      !
76      !
77      SELECT CASE( cd_type )        ! point of calculation
78      !
79      CASE( 'DYN' )                     ! T- and F-points
80         DO jk = jpkm1, 1, -1                ! pah1 at T-point
81            pah1(:,:,jk) = pahs1(:,:) * (  zratio + zc * ( 1._wp + TANH( - ( gdept_0(:,:,jk) - zh ) * zw) )  )
82         END DO
83         DO jk = jpkm1, 1, -1                ! pah2 at F-point (zdep2 is an approximation in zps-coord.)
84            DO jj = 1, jpjm1
85               DO ji = 1, jpim1
86                  zdep2 = (  gdept_0(ji,jj+1,jk) + gdept_0(ji+1,jj+1,jk)   &
87                     &     + gdept_0(ji,jj  ,jk) + gdept_0(ji+1,jj  ,jk)  ) * r1_4
88                  pah2(ji,jj,jk) = pahs2(ji,jj) * (  zratio + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  )
89               END DO
90            END DO
91         END DO
92         CALL lbc_lnk( 'ldfc1d_c2d', pah2, 'F', 1. )   ! Lateral boundary conditions
93         !
94      CASE( 'TRA' )                     ! U- and V-points (zdep1 & 2 are an approximation in zps-coord.)
95         DO jk = jpkm1, 1, -1
96            DO jj = 1, jpjm1
97               DO ji = 1, jpim1
98                  zdep1 = (  gdept_0(ji,jj,jk) + gdept_0(ji+1,jj,jk)  ) * 0.5_wp
99                  zdep2 = (  gdept_0(ji,jj,jk) + gdept_0(ji,jj+1,jk)  ) * 0.5_wp
100                  pah1(ji,jj,jk) = pahs1(ji,jj) * (  zratio + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) )  )
101                  pah2(ji,jj,jk) = pahs2(ji,jj) * (  zratio + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  )
102               END DO
103            END DO
104         END DO
105         ! Lateral boundary conditions
106         CALL lbc_lnk_multi( 'ldfc1d_c2d', pah1, 'U', 1. , pah2, 'V', 1. )   
107         !
108      CASE DEFAULT                        ! error
109         CALL ctl_stop( 'ldf_c1d: ', cd_type, ' Unknown, i.e. /= DYN or TRA' )
110      END SELECT
111      !
112   END SUBROUTINE ldf_c1d
113
114
115   SUBROUTINE ldf_c2d( cd_type, pUfac, knn, pah1, pah2 )
116      !!----------------------------------------------------------------------
117      !!                  ***  ROUTINE ldf_c2d  ***
118      !!             
119      !! ** Purpose :   2D eddy diffusivity/viscosity coefficients
120      !!
121      !! ** Method  :   2D eddy diffusivity coefficients F( e1 , e2 )
122      !!       laplacian   operator :   ah proportional to the scale factor      [m2/s]
123      !!       bilaplacian operator :   ah proportional to the (scale factor)^3  [m4/s]
124      !!       In both cases, pah0 is the maximum value reached by the coefficient
125      !!       at the Equator in case of e1=ra*rad= ~111km, not over the whole domain.
126      !!
127      !!   cd_type = TRA      pah1, pah2 defined at U- and V-points
128      !!             DYN      pah1, pah2 defined at T- and F-points
129      !!----------------------------------------------------------------------
130      CHARACTER(len=3)          , INTENT(in   ) ::   cd_type      ! DYNamique or TRAcers
131      REAL(wp)                  , INTENT(in   ) ::   pUfac        ! =1/2*Uc LAPlacian BiLaPlacian
132      INTEGER                   , INTENT(in   ) ::   knn          ! characteristic velocity   [m/s]
133      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pah1, pah2   ! eddy coefficients         [m2/s or m4/s]
134      !
135      INTEGER  ::   ji, jj, jk   ! dummy loop indices
136      INTEGER  ::   inn          ! local integer
137      !!----------------------------------------------------------------------
138      !
139      IF(lwp) THEN
140         WRITE(numout,*)
141         WRITE(numout,*) '   ldf_c2d :   aht = Ufac * max(e1,e2)   with Ufac = ', pUfac, ' m/s'
142         IF(lflush) CALL FLUSH(numout)
143      ENDIF
144      !
145      !
146      SELECT CASE( cd_type )        !==  surface values  ==!  (chosen grid point function of DYN or TRA)
147      !
148      CASE( 'DYN' )                       ! T- and F-points
149         DO jj = 1, jpj
150            DO ji = 1, jpi 
151               pah1(ji,jj,1) = pUfac * MAX( e1t(ji,jj) , e2t(ji,jj) )**knn
152               pah2(ji,jj,1) = pUfac * MAX( e1f(ji,jj) , e2f(ji,jj) )**knn
153            END DO
154         END DO
155      CASE( 'TRA' )                       ! U- and V-points
156         DO jj = 1, jpj 
157            DO ji = 1, jpi 
158               pah1(ji,jj,1) = pUfac * MAX( e1u(ji,jj), e2u(ji,jj) )**knn
159               pah2(ji,jj,1) = pUfac * MAX( e1v(ji,jj), e2v(ji,jj) )**knn
160            END DO
161         END DO
162      CASE DEFAULT                        ! error
163         CALL ctl_stop( 'ldf_c2d: ', cd_type, ' Unknown, i.e. /= DYN or TRA' )
164      END SELECT
165      !                             !==  deeper values = surface one  ==!  (except jpk)
166      DO jk = 2, jpkm1
167         pah1(:,:,jk) = pah1(:,:,1)
168         pah2(:,:,jk) = pah2(:,:,1)
169      END DO
170      !
171   END SUBROUTINE ldf_c2d
172
173   !!======================================================================
174END MODULE ldfc1d_c2d
Note: See TracBrowser for help on using the repository browser.