source: trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 7881

Last change on this file since 7881 was 7753, checked in by mocavero, 4 years ago

Reverting trunk to remove OpenMP

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