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

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 9366

Last change on this file since 9366 was 9366, checked in by andmirek, 6 years ago

#2050 first version. Compiled OK in moci test suite

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