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/UKMO/dev_r10448_bdyvol/src/OCE/ZDF – NEMO

source: NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE/ZDF/zdfddm.F90 @ 10455

Last change on this file since 10455 was 10068, checked in by nicolasmartin, 6 years ago

First part of modifications to have a common default header : fix typos and SVN keywords properties

  • Property svn:keywords set to Id
File size: 8.9 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, 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 indexocean time step
71      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm   !  Kz on momentum    (w-points)
72      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avt   !  Kz on temperature (w-points)
73      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avs   !  Kz on salinity    (w-points)
74      !
75      INTEGER  ::   ji, jj , jk     ! dummy loop indices
76      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
77      REAL(wp) ::   zdt, zds
78      REAL(wp) ::   zinr, zrr       !   -      -
79      REAL(wp) ::   zavft, zavfs    !   -      -
80      REAL(wp) ::   zavdt, zavds    !   -      -
81      REAL(wp), DIMENSION(jpi,jpj) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
82      !!----------------------------------------------------------------------
83      !
84      !                                                ! ===============
85      DO jk = 2, jpkm1                                 ! Horizontal slab
86         !                                             ! ===============
87         ! Define the mask
88         ! ---------------
89!!gm  WORK to be done:   change the code from vector optimisation to scalar one.
90!!gm                     ==>>>  test in the loop instead of use of mask arrays
91!!gm                            and many acces in memory
92         
93         DO jj = 1, jpj                !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
94            DO ji = 1, jpi
95               zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
96!!gm please, use e3w_n below
97                  &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) ) 
98               !
99               zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  )  &
100                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
101               zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  )  &
102                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
103               !
104               zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )
105               zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) 
106               IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
107               zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau
108            END DO
109         END DO
110
111         DO jj = 1, jpj                !==  indicators  ==!
112            DO ji = 1, jpi
113               ! stability indicator: msks=1 if rn2>0; 0 elsewhere
114               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
115               ELSE                                       ;   zmsks(ji,jj) = 1._wp
116               ENDIF
117               ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere           
118               IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp
119               ELSE                                       ;   zmskf(ji,jj) = 1._wp
120               ENDIF
121               ! diffusive layering indicators:
122               !     ! mskdl1=1 if 0< R <1; 0 elsewhere
123               IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp
124               ELSE                                       ;   zmskd1(ji,jj) = 1._wp
125               ENDIF
126               !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere
127               IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp
128               ELSE                                       ;   zmskd2(ji,jj) = 1._wp
129               ENDIF
130               !   mskdl3=1 if 0.5< R <1; 0 elsewhere
131               IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
132               ELSE                                                   ;   zmskd3(ji,jj) = 1._wp
133               ENDIF
134            END DO
135         END DO
136         ! mask zmsk in order to have avt and avs masked
137         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk)
138
139
140         ! Update avt and avs
141         ! ------------------
142         ! Constant eddy coefficient: reset to the background value
143         DO jj = 1, jpj
144            DO ji = 1, jpi
145               zinr = 1._wp / zrau(ji,jj)
146               ! salt fingering
147               zrr = zrau(ji,jj) / rn_hsbfr
148               zrr = zrr * zrr
149               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
150               zavft = 0.7 * zavfs * zinr
151               ! diffusive layering
152               zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
153               zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   &
154                  &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
155               ! add to the eddy viscosity coef. previously computed
156               p_avs(ji,jj,jk) = p_avt(ji,jj,jk) + zavfs + zavds
157               p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zavft + zavdt
158               p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
159            END DO
160         END DO
161         !                                                ! ===============
162      END DO                                              !   End of slab
163      !                                                   ! ===============
164      !
165      IF(ln_ctl) THEN
166         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', kdim=jpk)
167      ENDIF
168      !
169   END SUBROUTINE zdf_ddm
170   
171   !!======================================================================
172END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.