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

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ZDF/zdfddm.F90 @ 12236

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

Branch 2019/dev_r11943_MERGE_2019. Merge in changes from 2019/fix_sn_cfctl_ticket2328. Fully SETTE tested
  • Property svn:keywords set to Id
File size: 9.0 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
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   zdf_ddm       ! called by step.F90
30
31   !! * Substitutions
32#  include "vectopt_loop_substitute.h90"
33   !!----------------------------------------------------------------------
34   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
35   !! $Id$
36   !! Software governed by the CeCILL license (see ./LICENSE)
37   !!----------------------------------------------------------------------
38CONTAINS
39
40   SUBROUTINE zdf_ddm( kt, Kmm, p_avm, p_avt, p_avs )
41      !!----------------------------------------------------------------------
42      !!                  ***  ROUTINE zdf_ddm  ***
43      !!                   
44      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
45      !!              effect of salt fingering and diffusive convection.
46      !!
47      !! ** Method  :   Diapycnal mixing is increased in case of double
48      !!      diffusive mixing (i.e. salt fingering and diffusive layering)
49      !!      following Merryfield et al. (1999). The rate of double diffusive
50      !!      mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]):
51      !!         * salt fingering (Schmitt 1981):
52      !!      for R > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (R/rn_hsbfr)^6 )
53      !!      for R > 1 and rn2 > 0 : zavfs = O
54      !!      otherwise                : zavft = 0.7 zavs / R
55      !!         * diffusive layering (Federov 1988):
56      !!      for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) )
57      !!      otherwise                   : zavdt = 0
58      !!      for .5 < R < 1 and N^2 > 0 : zavds = zavdt (1.885 R -0.85)
59      !!      for  0 < R <.5 and N^2 > 0 : zavds = zavdt 0.15 R     
60      !!      otherwise                     : zavds = 0
61      !!         * update the eddy diffusivity:
62      !!      avt = avt + zavft + zavdt
63      !!      avs = avs + zavfs + zavds
64      !!      avm is required to remain at least above avt and avs.
65      !!     
66      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
67      !!
68      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
69      !!----------------------------------------------------------------------
70      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
71      INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index
72      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm   !  Kz on momentum    (w-points)
73      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avt   !  Kz on temperature (w-points)
74      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avs   !  Kz on salinity    (w-points)
75      !
76      INTEGER  ::   ji, jj , jk     ! dummy loop indices
77      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
78      REAL(wp) ::   zdt, zds
79      REAL(wp) ::   zinr, zrr       !   -      -
80      REAL(wp) ::   zavft, zavfs    !   -      -
81      REAL(wp) ::   zavdt, zavds    !   -      -
82      REAL(wp), DIMENSION(jpi,jpj) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
83      !!----------------------------------------------------------------------
84      !
85      !                                                ! ===============
86      DO jk = 2, jpkm1                                 ! Horizontal slab
87         !                                             ! ===============
88         ! Define the mask
89         ! ---------------
90!!gm  WORK to be done:   change the code from vector optimisation to scalar one.
91!!gm                     ==>>>  test in the loop instead of use of mask arrays
92!!gm                            and many acces in memory
93         
94         DO jj = 1, jpj                !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
95            DO ji = 1, jpi
96               zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
97!!gm please, use e3w(:,:,:,Kmm) below
98                  &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
99               !
100               zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  )  &
101                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
102               zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  )  &
103                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
104               !
105               zdt = zaw * ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) )
106               zds = zbw * ( ts(ji,jj,jk-1,jp_sal,Kmm) - ts(ji,jj,jk,jp_sal,Kmm) ) 
107               IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
108               zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau
109            END DO
110         END DO
111
112         DO jj = 1, jpj                !==  indicators  ==!
113            DO ji = 1, jpi
114               ! stability indicator: msks=1 if rn2>0; 0 elsewhere
115               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
116               ELSE                                       ;   zmsks(ji,jj) = 1._wp
117               ENDIF
118               ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere           
119               IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp
120               ELSE                                       ;   zmskf(ji,jj) = 1._wp
121               ENDIF
122               ! diffusive layering indicators:
123               !     ! mskdl1=1 if 0< R <1; 0 elsewhere
124               IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp
125               ELSE                                       ;   zmskd1(ji,jj) = 1._wp
126               ENDIF
127               !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere
128               IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp
129               ELSE                                       ;   zmskd2(ji,jj) = 1._wp
130               ENDIF
131               !   mskdl3=1 if 0.5< R <1; 0 elsewhere
132               IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
133               ELSE                                                   ;   zmskd3(ji,jj) = 1._wp
134               ENDIF
135            END DO
136         END DO
137         ! mask zmsk in order to have avt and avs masked
138         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk)
139
140
141         ! Update avt and avs
142         ! ------------------
143         ! Constant eddy coefficient: reset to the background value
144         DO jj = 1, jpj
145            DO ji = 1, jpi
146               zinr = 1._wp / zrau(ji,jj)
147               ! salt fingering
148               zrr = zrau(ji,jj) / rn_hsbfr
149               zrr = zrr * zrr
150               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
151               zavft = 0.7 * zavfs * zinr
152               ! diffusive layering
153               zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
154               zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   &
155                  &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
156               ! add to the eddy viscosity coef. previously computed
157               p_avs(ji,jj,jk) = p_avt(ji,jj,jk) + zavfs + zavds
158               p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zavft + zavdt
159               p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
160            END DO
161         END DO
162         !                                                ! ===============
163      END DO                                              !   End of slab
164      !                                                   ! ===============
165      !
166      IF(sn_cfctl%l_prtctl) THEN
167         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', kdim=jpk)
168      ENDIF
169      !
170   END SUBROUTINE zdf_ddm
171   
172   !!======================================================================
173END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.