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.
zdfric.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ZDF – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ZDF/zdfric.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.
  • Property svn:keywords set to Id
File size: 12.1 KB
Line 
1MODULE zdfric
2   !!======================================================================
3   !!                       ***  MODULE  zdfric  ***
4   !! Ocean physics:  vertical mixing coefficient compute from the local
5   !!                 Richardson number dependent formulation
6   !!======================================================================
7   !! History :  OPA  !  1987-09  (P. Andrich)  Original code
8   !!            4.0  !  1991-11  (G. Madec)
9   !!            7.0  !  1996-01  (G. Madec)  complete rewriting of multitasking suppression of common work arrays
10   !!            8.0  !  1997-06  (G. Madec)  complete rewriting of zdfmix
11   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
12   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
13   !!            3.3.1!  2011-09  (P. Oddo) Mixed layer depth parameterization
14   !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only
15   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!   zdf_ric_init  : initialization, namelist read, & parameters control
19   !!   zdf_ric       : update momentum and tracer Kz from the Richardson number
20   !!   ric_rst       : read/write RIC restart in ocean restart file
21   !!----------------------------------------------------------------------
22   USE oce            ! ocean dynamics and tracers variables
23   USE dom_oce        ! ocean space and time domain variables
24   USE zdf_oce        ! vertical physics: variables
25   USE phycst         ! physical constants
26   USE sbc_oce,  ONLY :   taum
27   !
28   USE in_out_manager ! I/O manager
29   USE iom            ! I/O manager library
30   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
31
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   zdf_ric         ! called by zdfphy.F90
37   PUBLIC   ric_rst         ! called by zdfphy.F90
38   PUBLIC   zdf_ric_init    ! called by nemogcm.F90
39
40   !                        !!* Namelist namzdf_ric : Richardson number dependent Kz *
41   INTEGER  ::   nn_ric      ! coefficient of the parameterization
42   REAL(wp) ::   rn_avmri    ! maximum value of the vertical eddy viscosity
43   REAL(wp) ::   rn_alp      ! coefficient of the parameterization
44   REAL(wp) ::   rn_ekmfc    ! Ekman Factor Coeff
45   REAL(wp) ::   rn_mldmin   ! minimum mixed layer (ML) depth   
46   REAL(wp) ::   rn_mldmax   ! maximum mixed layer depth
47   REAL(wp) ::   rn_wtmix    ! Vertical eddy Diff. in the ML
48   REAL(wp) ::   rn_wvmix    ! Vertical eddy Visc. in the ML
49   LOGICAL  ::   ln_mldw     ! Use or not the MLD parameters
50
51   !! * Substitutions
52#  include "vectopt_loop_substitute.h90"
53#  include "do_loop_substitute.h90"
54   !!----------------------------------------------------------------------
55   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
56   !! $Id$
57   !! Software governed by the CeCILL license (see ./LICENSE)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   SUBROUTINE zdf_ric_init
62      !!----------------------------------------------------------------------
63      !!                 ***  ROUTINE zdf_ric_init  ***
64      !!                   
65      !! ** Purpose :   Initialization of the vertical eddy diffusivity and
66      !!      viscosity coef. for the Richardson number dependent formulation.
67      !!
68      !! ** Method  :   Read the namzdf_ric namelist and check the parameter values
69      !!
70      !! ** input   :   Namelist namzdf_ric
71      !!
72      !! ** Action  :   increase by 1 the nstop flag is setting problem encounter
73      !!----------------------------------------------------------------------
74      INTEGER ::   ji, jj, jk   ! dummy loop indices
75      INTEGER ::   ios          ! Local integer output status for namelist read
76      !!
77      NAMELIST/namzdf_ric/ rn_avmri, rn_alp   , nn_ric  , rn_ekmfc,  &
78         &                rn_mldmin, rn_mldmax, rn_wtmix, rn_wvmix, ln_mldw
79      !!----------------------------------------------------------------------
80      !
81      READ  ( numnam_ref, namzdf_ric, IOSTAT = ios, ERR = 901)
82901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ric in reference namelist' )
83
84      READ  ( numnam_cfg, namzdf_ric, IOSTAT = ios, ERR = 902 )
85902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_ric in configuration namelist' )
86      IF(lwm) WRITE ( numond, namzdf_ric )
87      !
88      IF(lwp) THEN                   ! Control print
89         WRITE(numout,*)
90         WRITE(numout,*) 'zdf_ric_init : Ri depend vertical mixing scheme'
91         WRITE(numout,*) '~~~~~~~~~~~~'
92         WRITE(numout,*) '   Namelist namzdf_ric : set Kz=F(Ri) parameters'
93         WRITE(numout,*) '      maximum vertical viscosity        rn_avmri  = ', rn_avmri
94         WRITE(numout,*) '      coefficient                       rn_alp    = ', rn_alp
95         WRITE(numout,*) '      exponent                          nn_ric    = ', nn_ric
96         WRITE(numout,*) '      Ekman layer enhanced mixing       ln_mldw   = ', ln_mldw
97         WRITE(numout,*) '         Ekman Factor Coeff             rn_ekmfc  = ', rn_ekmfc
98         WRITE(numout,*) '         minimum mixed layer depth      rn_mldmin = ', rn_mldmin
99         WRITE(numout,*) '         maximum mixed layer depth      rn_mldmax = ', rn_mldmax
100         WRITE(numout,*) '         Vertical eddy Diff. in the ML  rn_wtmix  = ', rn_wtmix
101         WRITE(numout,*) '         Vertical eddy Visc. in the ML  rn_wvmix  = ', rn_wvmix
102      ENDIF
103      !
104      CALL ric_rst( nit000, 'READ' )  !* read or initialize all required files
105      !
106      IF( lwxios ) THEN
107         CALL iom_set_rstw_var_active('avt_k')
108         CALL iom_set_rstw_var_active('avm_k')
109      ENDIF
110   END SUBROUTINE zdf_ric_init
111
112
113   SUBROUTINE zdf_ric( kt, Kmm, p_sh2, p_avm, p_avt )
114      !!----------------------------------------------------------------------
115      !!                 ***  ROUTINE zdfric  ***
116      !!                   
117      !! ** Purpose :   Compute the before eddy viscosity and diffusivity as
118      !!                a function of the local richardson number.
119      !!
120      !! ** Method  :   Local richardson number dependent formulation of the
121      !!                vertical eddy viscosity and diffusivity coefficients.
122      !!                The eddy coefficients are given by:
123      !!                    avm = avm0 + avmb
124      !!                    avt = avm0 / (1 + rn_alp*ri)
125      !!                with ri  = N^2 / dz(u)**2
126      !!                         = e3w**2 * rn2/[ mi( dk(uu(:,:,:,Kbb)) )+mj( dk(vv(:,:,:,Kbb)) ) ]
127      !!                    avm0= rn_avmri / (1 + rn_alp*Ri)**nn_ric
128      !!                where ri is the before local Richardson number,
129      !!                rn_avmri is the maximum value reaches by avm and avt
130      !!                and rn_alp, nn_ric are adjustable parameters.
131      !!                Typical values : rn_alp=5. and nn_ric=2.
132      !!
133      !!      As second step compute Ekman depth from wind stress forcing
134      !!      and apply namelist provided vertical coeff within this depth.
135      !!      The Ekman depth is:
136      !!              Ustar = SQRT(Taum/rho0)
137      !!              ekd= rn_ekmfc * Ustar / f0
138      !!      Large et al. (1994, eq.24) suggest rn_ekmfc=0.7; however, the derivation
139      !!      of the above equation indicates the value is somewhat arbitrary; therefore
140      !!      we allow the freedom to increase or decrease this value, if the
141      !!      Ekman depth estimate appears too shallow or too deep, respectively.
142      !!      Ekd is then limited by rn_mldmin and rn_mldmax provided in the
143      !!      namelist
144      !!        N.B. the mask are required for implicit scheme, and surface
145      !!      and bottom value already set in zdfphy.F90
146      !!
147      !! ** Action  :   avm, avt  mixing coeff (inner domain values only)
148      !!
149      !! References : Pacanowski & Philander 1981, JPO, 1441-1451.
150      !!              PFJ Lermusiaux 2001.
151      !!----------------------------------------------------------------------
152      INTEGER                   , INTENT(in   ) ::   kt             ! ocean time-step
153      INTEGER                   , INTENT(in   ) ::   Kmm            ! ocean time level index
154      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term
155      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   ! momentum and tracer Kz (w-points)
156      !!
157      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
158      REAL(wp) ::   zcfRi, zav, zustar, zhek    ! local scalars
159      REAL(wp), DIMENSION(jpi,jpj) ::   zh_ekm  ! 2D workspace
160      !!----------------------------------------------------------------------
161      !
162      !                       !==  avm and avt = F(Richardson number)  ==!
163      DO_3D_10_10( 2, jpkm1 )
164         zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( p_sh2(ji,jj,jk) + 1.e-20 ) )  )
165         zav   = rn_avmri * zcfRi**nn_ric
166         !                          ! avm and avt coefficients
167         p_avm(ji,jj,jk) = MAX(  zav         , avmb(jk)  ) * wmask(ji,jj,jk)
168         p_avt(ji,jj,jk) = MAX(  zav * zcfRi , avtb(jk)  ) * wmask(ji,jj,jk)
169      END_3D
170      !
171!!gm BUG <<<<====  This param can't work at low latitude
172!!gm               it provides there much to thick mixed layer ( summer 150m in GYRE configuration !!! )
173      !
174      IF( ln_mldw ) THEN      !==  set a minimum value in the Ekman layer  ==!
175         !
176         DO_2D_00_00
177            zustar = SQRT( taum(ji,jj) * r1_rau0 )
178            zhek   = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )   ! Ekman depth
179            zh_ekm(ji,jj) = MAX(  rn_mldmin , MIN( zhek , rn_mldmax )  )   ! set allowed range
180         END_2D
181         DO_3D_00_00( 2, jpkm1 )
182            IF( gdept(ji,jj,jk,Kmm) < zh_ekm(ji,jj) ) THEN
183               p_avm(ji,jj,jk) = MAX(  p_avm(ji,jj,jk), rn_wvmix  ) * wmask(ji,jj,jk)
184               p_avt(ji,jj,jk) = MAX(  p_avt(ji,jj,jk), rn_wtmix  ) * wmask(ji,jj,jk)
185            ENDIF
186         END_3D
187      ENDIF
188      !
189   END SUBROUTINE zdf_ric
190
191
192   SUBROUTINE ric_rst( kt, cdrw )
193      !!---------------------------------------------------------------------
194      !!                   ***  ROUTINE ric_rst  ***
195      !!                     
196      !! ** Purpose :   Read or write TKE file (en) in restart file
197      !!
198      !! ** Method  :   use of IOM library
199      !!                if the restart does not contain TKE, en is either
200      !!                set to rn_emin or recomputed
201      !!----------------------------------------------------------------------
202      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
203      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
204      !
205      INTEGER ::   jit, jk    ! dummy loop indices
206      INTEGER ::   id1, id2   ! local integers
207      !!----------------------------------------------------------------------
208      !
209      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
210         !                                   ! ---------------
211         !           !* Read the restart file
212         IF( ln_rstart ) THEN
213            id1 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
214            id2 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
215            !
216            IF( MIN( id1, id2 ) > 0 ) THEN         ! restart exists => read it
217               CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k, ldxios = lrxios )
218               CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k, ldxios = lrxios )
219            ENDIF
220         ENDIF
221         !           !* otherwise Kz already set to the background value in zdf_phy_init
222         !
223      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
224         !                                   ! -------------------
225         IF(lwp) WRITE(numout,*) '---- ric-rst ----'
226         IF( lwxios ) CALL iom_swap(      cwxios_context          )
227         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k, ldxios = lwxios )
228         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k, ldxios = lwxios)
229         IF( lwxios ) CALL iom_swap(      cxios_context          )
230         !
231      ENDIF
232      !
233   END SUBROUTINE ric_rst
234
235   !!======================================================================
236END MODULE zdfric
Note: See TracBrowser for help on using the repository browser.