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/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 5951

Last change on this file since 5951 was 5951, checked in by timgraham, 8 years ago

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