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

source: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 3610

Last change on this file since 3610 was 3610, checked in by acc, 11 years ago

Branch dev_NOC_2012_r3555. #1006. Step 5: Merge in trunk changes between revision 3337 and 3385

  • Property svn:keywords set to Id
File size: 12.3 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   !!----------------------------------------------------------------------
10#if defined key_zdfddm   ||   defined key_esopa
11   !!----------------------------------------------------------------------
12   !!   'key_zdfddm' :                                     double diffusion
13   !!----------------------------------------------------------------------
14   !!   zdf_ddm       : compute the Ks for salinity
15   !!   zdf_ddm_init  : read namelist and control the parameters
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and tracers variables
18   USE dom_oce         ! ocean space and time domain variables
19   USE zdf_oce         ! ocean vertical physics variables
20   USE in_out_manager  ! I/O manager
21   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
22   USE prtctl          ! Print control
23   USE lib_mpp         ! MPP library
24   USE wrk_nemo        ! work arrays
25   USE timing          ! Timing
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   zdf_ddm       ! called by step.F90
31   PUBLIC   zdf_ddm_init  ! called by opa.F90
32   PUBLIC   zdf_ddm_alloc ! called by nemogcm.F90
33
34   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfddm = .TRUE.  !: double diffusive mixing flag
35
36   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avs    !: salinity vertical diffusivity coeff. at w-point
37   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   rrau   !: heat/salt buoyancy flux ratio
38
39   !                                  !!* Namelist namzdf_ddm : double diffusive mixing *
40   REAL(wp) ::   rn_avts  = 1.e-4_wp   ! maximum value of avs for salt fingering
41   REAL(wp) ::   rn_hsbfr = 1.6_wp     ! heat/salt buoyancy flux ratio
42
43   !! * Substitutions
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   INTEGER FUNCTION zdf_ddm_alloc()
53      !!----------------------------------------------------------------------
54      !!                ***  ROUTINE zdf_ddm_alloc  ***
55      !!----------------------------------------------------------------------
56      ALLOCATE( avs(jpi,jpj,jpk), rrau(jpi,jpj,jpk), STAT= zdf_ddm_alloc )
57      !
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: Rrau=alpha/beta dk[T]/dk[S]
74      !!      which is computed in rn2.F
75      !!         * salt fingering (Schmitt 1981):
76      !!      for Rrau > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (Rrau/rn_hsbfr)^6 )
77      !!      for Rrau > 1 and rn2 > 0 : zavfs = O
78      !!      otherwise                : zavft = 0.7 zavs / Rrau
79      !!         * diffusive layering (Federov 1988):
80      !!      for 0< Rrau < 1 and rn2 > 0 : zavdt = 1.3635e-6 
81      !!                                 * exp( 4.6 exp(-0.54 (1/Rrau-1) ) )
82      !!      otherwise                   : zavdt = 0
83      !!      for .5 < Rrau < 1 and rn2 > 0 : zavds = zavdt (1.885 Rrau -0.85)
84      !!      for  0 < Rrau <.5 and rn2 > 0 : zavds = zavdt 0.15 Rrau     
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) ::   zinr, zrr       ! temporary scalars
99      REAL(wp) ::   zavft, zavfs    !    -         -
100      REAL(wp) ::   zavdt, zavds    !    -         -
101      REAL(wp), POINTER, DIMENSION(:,:) ::   zmsks, zmskf, zmskd1, zmskd2, zmskd3
102      !!----------------------------------------------------------------------
103      !
104      IF( nn_timing == 1 )  CALL timing_start('zdf_ddm')
105      !
106      CALL wrk_alloc( jpi,jpj, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
107
108      !                                                ! ===============
109      DO jk = 2, jpkm1                                 ! Horizontal slab
110         !                                             ! ===============
111         ! Define the mask
112         ! ---------------
113         rrau(:,:,jk) = MAX( 1.e-20, rrau(:,:,jk) )         ! only retains positive value of rrau
114
115         DO jj = 1, jpj                                     ! indicators:
116            DO ji = 1, jpi
117               ! stability indicator: msks=1 if rn2>0; 0 elsewhere
118               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
119               ELSE                                       ;   zmsks(ji,jj) = 1._wp
120               ENDIF
121               ! salt fingering indicator: msksf=1 if rrau>1; 0 elsewhere           
122               IF( rrau(ji,jj,jk) <= 1.          ) THEN   ;   zmskf(ji,jj) = 0._wp
123               ELSE                                       ;   zmskf(ji,jj) = 1._wp
124               ENDIF
125               ! diffusive layering indicators:
126               !     ! mskdl1=1 if 0<rrau<1; 0 elsewhere
127               IF( rrau(ji,jj,jk) >= 1.          ) THEN   ;   zmskd1(ji,jj) = 0._wp
128               ELSE                                       ;   zmskd1(ji,jj) = 1._wp
129               ENDIF
130               !     ! mskdl2=1 if 0<rrau<0.5; 0 elsewhere
131               IF( rrau(ji,jj,jk) >= 0.5         ) THEN   ;   zmskd2(ji,jj) = 0._wp
132               ELSE                                       ;   zmskd2(ji,jj) = 1._wp
133               ENDIF
134               !   mskdl3=1 if 0.5<rrau<1; 0 elsewhere
135               IF( rrau(ji,jj,jk) <= 0.5 .OR. rrau(ji,jj,jk) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
136               ELSE                                                         ;   zmskd3(ji,jj) = 1._wp
137               ENDIF
138            END DO
139         END DO
140         ! mask zmsk in order to have avt and avs masked
141         zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk)
142
143
144         ! Update avt and avs
145         ! ------------------
146         ! Constant eddy coefficient: reset to the background value
147!CDIR NOVERRCHK
148         DO jj = 1, jpj
149!CDIR NOVERRCHK
150            DO ji = 1, jpi
151               zinr = 1./rrau(ji,jj,jk)
152               ! salt fingering
153               zrr = rrau(ji,jj,jk)/rn_hsbfr
154               zrr = zrr * zrr
155               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
156               zavft = 0.7 * zavfs * zinr
157               ! diffusive layering
158               zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
159               zavds = zavdt * zmsks(ji,jj) * (  (1.85 * rrau(ji,jj,jk) - 0.85 ) * zmskd3(ji,jj)   &
160                  &                            +  0.15 * rrau(ji,jj,jk)          * zmskd2(ji,jj)  )
161               ! add to the eddy viscosity coef. previously computed
162               avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds
163               avt (ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt
164               avm (ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
165            END DO
166         END DO
167
168
169         ! Increase avmu, avmv if necessary
170         ! --------------------------------
171!!gm to be changed following the definition of avm.
172         DO jj = 1, jpjm1
173            DO ji = 1, fs_jpim1   ! vector opt.
174               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    &
175                  &                  avt(ji,jj,jk), avt(ji+1,jj,jk),   &
176                  &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * umask(ji,jj,jk)
177               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    &
178                  &                  avt(ji,jj,jk), avt(ji,jj+1,jk),   &
179                  &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * vmask(ji,jj,jk)
180            END DO
181         END DO
182         !                                                ! ===============
183      END DO                                              !   End of slab
184      !                                                   ! ===============
185      !
186      CALL lbc_lnk( avt , 'W', 1._wp )     ! Lateral boundary conditions   (unchanged sign)
187      CALL lbc_lnk( avs , 'W', 1._wp )
188      CALL lbc_lnk( avm , 'W', 1._wp )
189      CALL lbc_lnk( avmu, 'U', 1._wp ) 
190      CALL lbc_lnk( avmv, 'V', 1._wp )
191
192      IF(ln_ctl) THEN
193         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
194         CALL prt_ctl(tab3d_1=avmu, clinfo1=' ddm  - u: ', mask1=umask, &
195            &         tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk)
196      ENDIF
197      !
198      CALL wrk_dealloc( jpi,jpj, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
199      !
200      IF( nn_timing == 1 )  CALL timing_stop('zdf_ddm')
201      !
202   END SUBROUTINE zdf_ddm
203   
204   
205   SUBROUTINE zdf_ddm_init
206      !!----------------------------------------------------------------------
207      !!                  ***  ROUTINE zdf_ddm_init  ***
208      !!
209      !! ** Purpose :   Initialization of double diffusion mixing scheme
210      !!
211      !! ** Method  :   Read the namzdf_ddm namelist and check the parameter values
212      !!              called by zdf_ddm at the first timestep (nit000)
213      !!----------------------------------------------------------------------
214      NAMELIST/namzdf_ddm/ rn_avts, rn_hsbfr
215      !!----------------------------------------------------------------------
216      !
217      REWIND( numnam )                ! Read Namelist namzdf_ddm : double diffusion mixing scheme
218      READ  ( numnam, namzdf_ddm )
219      !
220      IF(lwp) THEN                    ! Parameter print
221         WRITE(numout,*)
222         WRITE(numout,*) 'zdf_ddm : double diffusive mixing'
223         WRITE(numout,*) '~~~~~~~'
224         WRITE(numout,*) '   Namelist namzdf_ddm : set dd mixing parameter'
225         WRITE(numout,*) '      maximum avs for dd mixing      rn_avts   = ', rn_avts
226         WRITE(numout,*) '      heat/salt buoyancy flux ratio  rn_hsbfr  = ', rn_hsbfr
227      ENDIF
228      !
229      !                               ! allocate zdfddm arrays
230      IF( zdf_ddm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ddm_init : unable to allocate arrays' )
231      !                               ! initialization to masked Kz
232      avs(:,:,:) = rn_avt0 * tmask(:,:,:) 
233      !
234   END SUBROUTINE zdf_ddm_init
235
236#else
237   !!----------------------------------------------------------------------
238   !!   Default option :          Dummy module          No double diffusion
239   !!----------------------------------------------------------------------
240   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfddm = .FALSE.   !: double diffusion flag
241CONTAINS
242   SUBROUTINE zdf_ddm( kt )           ! Dummy routine
243      WRITE(*,*) 'zdf_ddm: You should not have seen this print! error?', kt
244   END SUBROUTINE zdf_ddm
245   SUBROUTINE zdf_ddm_init            ! Dummy routine
246   END SUBROUTINE zdf_ddm_init
247#endif
248
249   !!======================================================================
250END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.