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/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 9019

Last change on this file since 9019 was 9019, checked in by timgraham, 6 years ago

Merge of dev_CNRS_2017 into branch

  • 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   USE timing         ! Timing
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   zdf_ddm       ! called by step.F90
31
32   !! * Substitutions
33#  include "vectopt_loop_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( 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   ) ::   kt       ! ocean time-step indexocean time step
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      IF( ln_timing )   CALL timing_start('zdf_ddm')
86      !
87      !                                                ! ===============
88      DO jk = 2, jpkm1                                 ! Horizontal slab
89         !                                             ! ===============
90         ! Define the mask
91         ! ---------------
92!!gm  WORK to be done:   change the code from vector optimisation to scalar one.
93!!gm                     ==>>>  test in the loop instead of use of mask arrays
94!!gm                            and many acces in memory
95         
96         DO jj = 1, jpj                !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
97            DO ji = 1, jpi
98               zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
99!!gm please, use e3w_n below
100                  &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) ) 
101               !
102               zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  )  &
103                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
104               zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  )  &
105                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
106               !
107               zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )
108               zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) 
109               IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
110               zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau
111            END DO
112         END DO
113
114         DO jj = 1, jpj                !==  indicators  ==!
115            DO ji = 1, jpi
116               ! stability indicator: msks=1 if rn2>0; 0 elsewhere
117               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
118               ELSE                                       ;   zmsks(ji,jj) = 1._wp
119               ENDIF
120               ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere           
121               IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp
122               ELSE                                       ;   zmskf(ji,jj) = 1._wp
123               ENDIF
124               ! diffusive layering indicators:
125               !     ! mskdl1=1 if 0< R <1; 0 elsewhere
126               IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp
127               ELSE                                       ;   zmskd1(ji,jj) = 1._wp
128               ENDIF
129               !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere
130               IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp
131               ELSE                                       ;   zmskd2(ji,jj) = 1._wp
132               ENDIF
133               !   mskdl3=1 if 0.5< R <1; 0 elsewhere
134               IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
135               ELSE                                                   ;   zmskd3(ji,jj) = 1._wp
136               ENDIF
137            END DO
138         END DO
139         ! mask zmsk in order to have avt and avs masked
140         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk)
141
142
143         ! Update avt and avs
144         ! ------------------
145         ! Constant eddy coefficient: reset to the background value
146         DO jj = 1, jpj
147            DO ji = 1, jpi
148               zinr = 1._wp / zrau(ji,jj)
149               ! salt fingering
150               zrr = zrau(ji,jj) / rn_hsbfr
151               zrr = zrr * zrr
152               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
153               zavft = 0.7 * zavfs * zinr
154               ! diffusive layering
155               zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
156               zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   &
157                  &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
158               ! add to the eddy viscosity coef. previously computed
159               p_avs(ji,jj,jk) = p_avt(ji,jj,jk) + zavfs + zavds
160               p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zavft + zavdt
161               p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
162            END DO
163         END DO
164         !                                                ! ===============
165      END DO                                              !   End of slab
166      !                                                   ! ===============
167      !
168      IF(ln_ctl) THEN
169         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
170      ENDIF
171      !
172      IF( ln_timing )   CALL timing_stop('zdf_ddm')
173      !
174   END SUBROUTINE zdf_ddm
175   
176   !!======================================================================
177END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.