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

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago

All wrk_alloc removed

  • Property svn:keywords set to Id
File size: 13.5 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   !!----------------------------------------------------------------------
11#if defined key_zdfddm
12   !!----------------------------------------------------------------------
13   !!   'key_zdfddm' :                                     double diffusion
14   !!----------------------------------------------------------------------
15   !!   zdf_ddm       : compute the Ks for salinity
16   !!   zdf_ddm_init  : read namelist and control the parameters
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers variables
19   USE dom_oce         ! ocean space and time domain variables
20   USE zdf_oce         ! ocean vertical physics variables
21   USE eosbn2         ! equation of state
22   !
23   USE in_out_manager  ! I/O manager
24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
25   USE prtctl          ! Print control
26   USE lib_mpp         ! MPP library
27   USE timing          ! Timing
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   zdf_ddm       ! called by step.F90
33   PUBLIC   zdf_ddm_init  ! called by opa.F90
34   PUBLIC   zdf_ddm_alloc ! called by nemogcm.F90
35
36   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfddm = .TRUE.  !: double diffusive mixing flag
37
38   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avs   !: salinity vertical diffusivity coeff. at w-point
39
40   !                       !!* Namelist namzdf_ddm : double diffusive mixing *
41   REAL(wp) ::   rn_avts    ! maximum value of avs for salt fingering
42   REAL(wp) ::   rn_hsbfr   ! heat/salt buoyancy flux ratio
43
44   !! * Substitutions
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
48   !! $Id$
49   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53   INTEGER FUNCTION zdf_ddm_alloc()
54      !!----------------------------------------------------------------------
55      !!                ***  ROUTINE zdf_ddm_alloc  ***
56      !!----------------------------------------------------------------------
57      ALLOCATE( avs(jpi,jpj,jpk) , STAT= zdf_ddm_alloc )
58      IF( lk_mpp             )   CALL mpp_sum ( zdf_ddm_alloc )
59      IF( zdf_ddm_alloc /= 0 )   CALL ctl_warn('zdf_ddm_alloc: failed to allocate arrays')
60   END FUNCTION zdf_ddm_alloc
61
62
63   SUBROUTINE zdf_ddm( kt )
64      !!----------------------------------------------------------------------
65      !!                  ***  ROUTINE zdf_ddm  ***
66      !!                   
67      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
68      !!              effect of salt fingering and diffusive convection.
69      !!
70      !! ** Method  :   Diapycnal mixing is increased in case of double
71      !!      diffusive mixing (i.e. salt fingering and diffusive layering)
72      !!      following Merryfield et al. (1999). The rate of double diffusive
73      !!      mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]):
74      !!         * salt fingering (Schmitt 1981):
75      !!      for R > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (R/rn_hsbfr)^6 )
76      !!      for R > 1 and rn2 > 0 : zavfs = O
77      !!      otherwise                : zavft = 0.7 zavs / R
78      !!         * diffusive layering (Federov 1988):
79      !!      for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) )
80      !!      otherwise                   : zavdt = 0
81      !!      for .5 < R < 1 and N^2 > 0 : zavds = zavdt (1.885 R -0.85)
82      !!      for  0 < R <.5 and N^2 > 0 : zavds = zavdt 0.15 R     
83      !!      otherwise                     : zavds = 0
84      !!         * update the eddy diffusivity:
85      !!      avt = avt + zavft + zavdt
86      !!      avs = avs + zavfs + zavds
87      !!      avmu, avmv are required to remain at least above avt and avs.
88      !!     
89      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
90      !!
91      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
92      !!----------------------------------------------------------------------
93      INTEGER, INTENT(in) ::   kt   ! ocean time-step indexocean time step
94      !
95      INTEGER  ::   ji, jj , jk     ! dummy loop indices
96      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
97      REAL(wp) ::   zdt, zds
98      REAL(wp) ::   zinr, zrr       !   -      -
99      REAL(wp) ::   zavft, zavfs    !   -      -
100      REAL(wp) ::   zavdt, zavds    !   -      -
101      REAL(wp), DIMENSION(jpi,jpj) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
102      !!----------------------------------------------------------------------
103      !
104      IF( nn_timing == 1 )  CALL timing_start('zdf_ddm')
105      !
106      !
107      !                                                ! ===============
108      DO jk = 2, jpkm1                                 ! Horizontal slab
109         !                                             ! ===============
110         ! Define the mask
111         ! ---------------
112         DO jj = 1, jpj                                ! R=zrau = (alpha / beta) (dk[t] / dk[s])
113            DO ji = 1, jpi
114               zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
115                  &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) ) 
116               !
117               zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  )  &
118                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
119               zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  )  &
120                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
121               !
122               zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )
123               zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) 
124               IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
125               zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau
126            END DO
127         END DO
128
129         DO jj = 1, jpj                                     ! indicators:
130            DO ji = 1, jpi
131               ! stability indicator: msks=1 if rn2>0; 0 elsewhere
132               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
133               ELSE                                       ;   zmsks(ji,jj) = 1._wp
134               ENDIF
135               ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere           
136               IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp
137               ELSE                                       ;   zmskf(ji,jj) = 1._wp
138               ENDIF
139               ! diffusive layering indicators:
140               !     ! mskdl1=1 if 0< R <1; 0 elsewhere
141               IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp
142               ELSE                                       ;   zmskd1(ji,jj) = 1._wp
143               ENDIF
144               !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere
145               IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp
146               ELSE                                       ;   zmskd2(ji,jj) = 1._wp
147               ENDIF
148               !   mskdl3=1 if 0.5< R <1; 0 elsewhere
149               IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
150               ELSE                                                   ;   zmskd3(ji,jj) = 1._wp
151               ENDIF
152            END DO
153         END DO
154         ! mask zmsk in order to have avt and avs masked
155         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk)
156
157
158         ! Update avt and avs
159         ! ------------------
160         ! Constant eddy coefficient: reset to the background value
161         DO jj = 1, jpj
162            DO ji = 1, jpi
163               zinr = 1._wp / zrau(ji,jj)
164               ! salt fingering
165               zrr = zrau(ji,jj) / rn_hsbfr
166               zrr = zrr * zrr
167               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
168               zavft = 0.7 * zavfs * zinr
169               ! diffusive layering
170               zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
171               zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   &
172                  &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
173               ! add to the eddy viscosity coef. previously computed
174# if defined key_zdftmx_new
175               ! key_zdftmx_new: New internal wave-driven param: use avs value computed by zdftmx
176               avs (ji,jj,jk) = avs(ji,jj,jk) + zavfs + zavds
177# else
178               avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds
179# endif
180               avt (ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt
181               avm (ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
182            END DO
183         END DO
184
185
186         ! Increase avmu, avmv if necessary
187         ! --------------------------------
188!!gm to be changed following the definition of avm.
189         DO jj = 1, jpjm1
190            DO ji = 1, fs_jpim1   ! vector opt.
191               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    &
192                  &                  avt(ji,jj,jk), avt(ji+1,jj,jk),   &
193                  &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * wumask(ji,jj,jk)
194               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    &
195                  &                  avt(ji,jj,jk), avt(ji,jj+1,jk),   &
196                  &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * wvmask(ji,jj,jk)
197            END DO
198         END DO
199         !                                                ! ===============
200      END DO                                              !   End of slab
201      !                                                   ! ===============
202      !
203      CALL lbc_lnk( avt , 'W', 1._wp )     ! Lateral boundary conditions   (unchanged sign)
204      CALL lbc_lnk( avs , 'W', 1._wp )
205      CALL lbc_lnk( avm , 'W', 1._wp )
206      CALL lbc_lnk( avmu, 'U', 1._wp ) 
207      CALL lbc_lnk( avmv, 'V', 1._wp )
208
209      IF(ln_ctl) THEN
210         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
211         CALL prt_ctl(tab3d_1=avmu, clinfo1=' ddm  - u: ', mask1=umask, &
212            &         tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk)
213      ENDIF
214      !
215      !
216      IF( nn_timing == 1 )  CALL timing_stop('zdf_ddm')
217      !
218   END SUBROUTINE zdf_ddm
219   
220   
221   SUBROUTINE zdf_ddm_init
222      !!----------------------------------------------------------------------
223      !!                  ***  ROUTINE zdf_ddm_init  ***
224      !!
225      !! ** Purpose :   Initialization of double diffusion mixing scheme
226      !!
227      !! ** Method  :   Read the namzdf_ddm namelist and check the parameter values
228      !!              called by zdf_ddm at the first timestep (nit000)
229      !!----------------------------------------------------------------------
230      INTEGER ::   ios   ! local integer
231      !!
232      NAMELIST/namzdf_ddm/ rn_avts, rn_hsbfr
233      !!----------------------------------------------------------------------
234      !
235      REWIND( numnam_ref )              ! Namelist namzdf_ddm in reference namelist : Double diffusion mixing scheme
236      READ  ( numnam_ref, namzdf_ddm, IOSTAT = ios, ERR = 901)
237901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ddm in reference namelist', lwp )
238
239      REWIND( numnam_cfg )              ! Namelist namzdf_ddm in configuration namelist : Double diffusion mixing scheme
240      READ  ( numnam_cfg, namzdf_ddm, IOSTAT = ios, ERR = 902 )
241902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ddm in configuration namelist', lwp )
242      IF(lwm) WRITE ( numond, namzdf_ddm )
243      !
244      IF(lwp) THEN                    ! Parameter print
245         WRITE(numout,*)
246         WRITE(numout,*) 'zdf_ddm : double diffusive mixing'
247         WRITE(numout,*) '~~~~~~~'
248         WRITE(numout,*) '   Namelist namzdf_ddm : set dd mixing parameter'
249         WRITE(numout,*) '      maximum avs for dd mixing      rn_avts   = ', rn_avts
250         WRITE(numout,*) '      heat/salt buoyancy flux ratio  rn_hsbfr  = ', rn_hsbfr
251      ENDIF
252      !
253      !                               ! allocate zdfddm arrays
254      IF( zdf_ddm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ddm_init : unable to allocate arrays' )
255      !                               ! initialization to masked Kz
256      avs(:,:,:) = rn_avt0 * wmask(:,:,:) 
257      !
258   END SUBROUTINE zdf_ddm_init
259
260#else
261   !!----------------------------------------------------------------------
262   !!   Default option :          Dummy module          No double diffusion
263   !!----------------------------------------------------------------------
264   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfddm = .FALSE.   !: double diffusion flag
265CONTAINS
266   SUBROUTINE zdf_ddm( kt )           ! Dummy routine
267      WRITE(*,*) 'zdf_ddm: You should not have seen this print! error?', kt
268   END SUBROUTINE zdf_ddm
269   SUBROUTINE zdf_ddm_init            ! Dummy routine
270   END SUBROUTINE zdf_ddm_init
271#endif
272
273   !!======================================================================
274END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.