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.
zdfddm.F90 in branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 8055

Last change on this file since 8055 was 8055, checked in by gm, 7 years ago

wrk_OMP: 2nd step: add OMP processor distribution in ZDF package
  • Property svn:keywords set to Id
File size: 9.2 KB
Line 
1MODULE zdfddm
2   !!======================================================================
3   !!                       ***  MODULE  zdfddm  ***
4   !! Ocean physics : double diffusion mixing parameterization
5   !!======================================================================
6   !! History :  OPA  ! 2000-08  (G. Madec)  double diffusive mixing
7   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
8   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
9   !!            3.6  ! 2013-04  (G. Madec, F. Roquet) zrau compute locally using interpolation of alpha & beta
10   !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   zdf_ddm       : compute the Kz for salinity
15   !!----------------------------------------------------------------------
16   USE oce            ! ocean dynamics and tracers variables
17   USE dom_oce        ! ocean space and time domain variables
18   USE zdf_oce        ! ocean vertical physics variables
19   USE eosbn2         ! equation of state
20   !
21   USE in_out_manager ! I/O manager
22   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
23   USE prtctl         ! Print control
24   USE lib_mpp        ! MPP library
25   USE timing         ! Timing
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   zdf_ddm       ! called by step.F90
31
32   !! * Substitutions
33#  include "domain_substitute.h90"   
34   !!----------------------------------------------------------------------
35   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
36   !! $Id$
37   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE zdf_ddm( ARG_2D, kt, p_avm, p_avt, p_avs )
42      !!----------------------------------------------------------------------
43      !!                  ***  ROUTINE zdf_ddm  ***
44      !!                   
45      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
46      !!              effect of salt fingering and diffusive convection.
47      !!
48      !! ** Method  :   Diapycnal mixing is increased in case of double
49      !!      diffusive mixing (i.e. salt fingering and diffusive layering)
50      !!      following Merryfield et al. (1999). The rate of double diffusive
51      !!      mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]):
52      !!         * salt fingering (Schmitt 1981):
53      !!      for R > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (R/rn_hsbfr)^6 )
54      !!      for R > 1 and rn2 > 0 : zavfs = O
55      !!      otherwise                : zavft = 0.7 zavs / R
56      !!         * diffusive layering (Federov 1988):
57      !!      for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) )
58      !!      otherwise                   : zavdt = 0
59      !!      for .5 < R < 1 and N^2 > 0 : zavds = zavdt (1.885 R -0.85)
60      !!      for  0 < R <.5 and N^2 > 0 : zavds = zavdt 0.15 R     
61      !!      otherwise                     : zavds = 0
62      !!         * update the eddy diffusivity:
63      !!      avt = avt + zavft + zavdt
64      !!      avs = avs + zavfs + zavds
65      !!      avm is required to remain at least above avt and avs.
66      !!     
67      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
68      !!
69      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
70      !!----------------------------------------------------------------------
71      INTEGER, INTENT(in   ) ::   ARG_2D   ! inner domain start-end i-indices
72      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step indexocean time step
73      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm   !  Kz on momentum    (w-points)
74      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avt   !  Kz on temperature (w-points)
75      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avs   !  Kz on salinity    (w-points)
76      !
77      INTEGER  ::   ji, jj , jk     ! dummy loop indices
78      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
79      REAL(wp) ::   zdt, zds
80      REAL(wp) ::   zinr, zrr       !   -      -
81      REAL(wp) ::   zavft, zavfs    !   -      -
82      REAL(wp) ::   zavdt, zavds    !   -      -
83      REAL(wp), DIMENSION(WRK_2D) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
84      !!----------------------------------------------------------------------
85      !
86      IF( nn_timing == 1 )   CALL timing_start('zdf_ddm')
87      !
88      !                                                ! ===============
89      DO jk = 2, jpkm1                                 ! Horizontal slab
90         !                                             ! ===============
91         ! Define the mask
92         ! ---------------
93!!gm  WORK to be done:   change the code from vector optimisation to scalar one.
94!!gm                     ==>>>  test in the loop instead of use of mask arrays
95!!gm                            and many acces in memory
96         
97         DO jj = k_Jstr, k_Jend        !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
98            DO ji = k_Jstr, k_Iend
99               zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
100!!gm please, use e3w_n below
101                  &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) ) 
102               !
103               zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  )  &
104                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
105               zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  )  &
106                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
107               !
108               zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )
109               zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) 
110               IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
111               zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau
112            END DO
113         END DO
114
115         DO jj = k_Jstr, k_Jend          !==  indicators  ==!
116            DO ji = k_Jstr, k_Iend
117               ! stability indicator: msks=1 if rn2>0; 0 elsewhere
118               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
119               ELSE                                       ;   zmsks(ji,jj) = 1._wp
120               ENDIF
121               ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere           
122               IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp
123               ELSE                                       ;   zmskf(ji,jj) = 1._wp
124               ENDIF
125               ! diffusive layering indicators:
126               !     ! mskdl1=1 if 0< R <1; 0 elsewhere
127               IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp
128               ELSE                                       ;   zmskd1(ji,jj) = 1._wp
129               ENDIF
130               !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere
131               IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp
132               ELSE                                       ;   zmskd2(ji,jj) = 1._wp
133               ENDIF
134               !   mskdl3=1 if 0.5< R <1; 0 elsewhere
135               IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
136               ELSE                                                   ;   zmskd3(ji,jj) = 1._wp
137               ENDIF
138            END DO
139         END DO
140         ! mask zmsk in order to have avt and avs masked
141         zmsks(WRK_2D) = zmsks(WRK_2D) * wmask(WRK_2D,jk)
142
143
144         ! Update avt and avs
145         ! ------------------
146         ! Constant eddy coefficient: reset to the background value
147         DO jj = k_Jstr, k_Jend
148            DO ji = k_Jstr, k_Iend
149               zinr = 1._wp / zrau(ji,jj)
150               ! salt fingering
151               zrr = zrau(ji,jj) / rn_hsbfr
152               zrr = zrr * zrr
153               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
154               zavft = 0.7 * zavfs * zinr
155               ! diffusive layering
156               zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
157               zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   &
158                  &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
159               ! add to the eddy viscosity coef. previously computed
160               p_avs(ji,jj,jk) = p_avt(ji,jj,jk) + zavfs + zavds
161               p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zavft + zavdt
162               p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
163            END DO
164         END DO
165         !                                                ! ===============
166      END DO                                              !   End of slab
167      !                                                   ! ===============
168      !
169      IF(ln_ctl) THEN
170         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
171      ENDIF
172      !
173      IF( nn_timing == 1 )  CALL timing_stop('zdf_ddm')
174      !
175   END SUBROUTINE zdf_ddm
176   
177   !!======================================================================
178END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.