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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 12.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   !!----------------------------------------------------------------------
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
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   zdf_ddm       ! called by step.F90
29   PUBLIC   zdf_ddm_init  ! called by opa.F90
30   PUBLIC   zdf_ddm_alloc ! called by nemogcm.F90
31
32   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfddm = .TRUE.  !: double diffusive mixing flag
33
34   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avs    !: salinity vertical diffusivity coeff. at w-point
35   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   rrau   !: heat/salt buoyancy flux ratio
36
37   !                                  !!* Namelist namzdf_ddm : double diffusive mixing *
38   REAL(wp) ::   rn_avts  = 1.e-4_wp   ! maximum value of avs for salt fingering
39   REAL(wp) ::   rn_hsbfr = 1.6_wp     ! heat/salt buoyancy flux ratio
40
41   !! * Control permutation of array indices
42#  include "zdfddm_ftrans.h90"
43#  include "oce_ftrans.h90"
44#  include "dom_oce_ftrans.h90"
45#  include "zdf_oce_ftrans.h90"
46
47   !! * Substitutions
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
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), rrau(jpi,jpj,jpk), STAT= zdf_ddm_alloc )
61      !
62      IF( lk_mpp             )   CALL mpp_sum ( zdf_ddm_alloc )
63      IF( zdf_ddm_alloc /= 0 )   CALL ctl_warn('zdf_ddm_alloc: failed to allocate arrays')
64   END FUNCTION zdf_ddm_alloc
65
66
67   SUBROUTINE zdf_ddm( kt )
68      !!----------------------------------------------------------------------
69      !!                  ***  ROUTINE zdf_ddm  ***
70      !!                   
71      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
72      !!              effect of salt fingering and diffusive convection.
73      !!
74      !! ** Method  :   Diapycnal mixing is increased in case of double
75      !!      diffusive mixing (i.e. salt fingering and diffusive layering)
76      !!      following Merryfield et al. (1999). The rate of double diffusive
77      !!      mixing depend on the buoyancy ratio: Rrau=alpha/beta dk[T]/dk[S]
78      !!      which is computed in rn2.F
79      !!         * salt fingering (Schmitt 1981):
80      !!      for Rrau > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (Rrau/rn_hsbfr)^6 )
81      !!      for Rrau > 1 and rn2 > 0 : zavfs = O
82      !!      otherwise                : zavft = 0.7 zavs / Rrau
83      !!         * diffusive layering (Federov 1988):
84      !!      for 0< Rrau < 1 and rn2 > 0 : zavdt = 1.3635e-6 
85      !!                                 * exp( 4.6 exp(-0.54 (1/Rrau-1) ) )
86      !!      otherwise                   : zavdt = 0
87      !!      for .5 < Rrau < 1 and rn2 > 0 : zavds = zavdt (1.885 Rrau -0.85)
88      !!      for  0 < Rrau <.5 and rn2 > 0 : zavds = zavdt 0.15 Rrau     
89      !!      otherwise                     : zavds = 0
90      !!         * update the eddy diffusivity:
91      !!      avt = avt + zavft + zavdt
92      !!      avs = avs + zavfs + zavds
93      !!      avmu, avmv are required to remain at least above avt and avs.
94      !!     
95      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
96      !!
97      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
98      !!----------------------------------------------------------------------
99      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
100      USE wrk_nemo, ONLY:   zmsks  => wrk_2d_1 , zmskf  => wrk_2d_2 , zmskd1 => wrk_2d_3   ! 2D workspace
101      USE wrk_nemo, ONLY:   zmskd2 => wrk_2d_4 , zmskd3 => wrk_2d_5                        !  -      -
102      !
103      INTEGER, INTENT(in) ::   kt   ! ocean time-step indexocean time step
104      !
105      INTEGER  ::   ji, jj , jk     ! dummy loop indices
106      REAL(wp) ::   zinr, zrr       ! temporary scalars
107      REAL(wp) ::   zavft, zavfs    !    -         -
108      REAL(wp) ::   zavdt, zavds    !    -         -
109      !!----------------------------------------------------------------------
110
111      IF( wrk_in_use(2, 1,2,3,4,5) ) THEN
112         CALL ctl_stop('zdf_ddm: Requested workspace arrays already in use')   ;   RETURN
113      ENDIF
114
115      !                                                ! ===============
116      DO jk = 2, jpkm1                                 ! Horizontal slab
117         !                                             ! ===============
118         ! Define the mask
119         ! ---------------
120         rrau(:,:,jk) = MAX( 1.e-20, rrau(:,:,jk) )         ! only retains positive value of rrau
121
122         DO jj = 1, jpj                                     ! indicators:
123            DO ji = 1, jpi
124               ! stability indicator: msks=1 if rn2>0; 0 elsewhere
125               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
126               ELSE                                       ;   zmsks(ji,jj) = 1._wp
127               ENDIF
128               ! salt fingering indicator: msksf=1 if rrau>1; 0 elsewhere           
129               IF( rrau(ji,jj,jk) <= 1.          ) THEN   ;   zmskf(ji,jj) = 0._wp
130               ELSE                                       ;   zmskf(ji,jj) = 1._wp
131               ENDIF
132               ! diffusive layering indicators:
133               !     ! mskdl1=1 if 0<rrau<1; 0 elsewhere
134               IF( rrau(ji,jj,jk) >= 1.          ) THEN   ;   zmskd1(ji,jj) = 0._wp
135               ELSE                                       ;   zmskd1(ji,jj) = 1._wp
136               ENDIF
137               !     ! mskdl2=1 if 0<rrau<0.5; 0 elsewhere
138               IF( rrau(ji,jj,jk) >= 0.5         ) THEN   ;   zmskd2(ji,jj) = 0._wp
139               ELSE                                       ;   zmskd2(ji,jj) = 1._wp
140               ENDIF
141               !   mskdl3=1 if 0.5<rrau<1; 0 elsewhere
142               IF( rrau(ji,jj,jk) <= 0.5 .OR. rrau(ji,jj,jk) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
143               ELSE                                                         ;   zmskd3(ji,jj) = 1._wp
144               ENDIF
145            END DO
146         END DO
147         ! mask zmsk in order to have avt and avs masked
148         zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk)
149
150
151         ! Update avt and avs
152         ! ------------------
153         ! Constant eddy coefficient: reset to the background value
154!CDIR NOVERRCHK
155         DO jj = 1, jpj
156!CDIR NOVERRCHK
157            DO ji = 1, jpi
158               zinr = 1./rrau(ji,jj,jk)
159               ! salt fingering
160               zrr = rrau(ji,jj,jk)/rn_hsbfr
161               zrr = zrr * zrr
162               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
163               zavft = 0.7 * zavfs * zinr
164               ! diffusive layering
165               zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
166               zavds = zavdt * zmsks(ji,jj) * (  (1.85 * rrau(ji,jj,jk) - 0.85 ) * zmskd3(ji,jj)   &
167                  &                            +  0.15 * rrau(ji,jj,jk)          * zmskd2(ji,jj)  )
168               ! add to the eddy viscosity coef. previously computed
169               avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds
170               avt (ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt
171               avm (ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
172            END DO
173         END DO
174
175
176         ! Increase avmu, avmv if necessary
177         ! --------------------------------
178!!gm to be changed following the definition of avm.
179         DO jj = 1, jpjm1
180            DO ji = 1, fs_jpim1   ! vector opt.
181               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    &
182                  &                  avt(ji,jj,jk), avt(ji+1,jj,jk),   &
183                  &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * umask(ji,jj,jk)
184               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    &
185                  &                  avt(ji,jj,jk), avt(ji,jj+1,jk),   &
186                  &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * vmask(ji,jj,jk)
187            END DO
188         END DO
189         !                                                ! ===============
190      END DO                                              !   End of slab
191      !                                                   ! ===============
192      !
193      CALL lbc_lnk( avt , 'W', 1._wp )     ! Lateral boundary conditions   (unchanged sign)
194      CALL lbc_lnk( avs , 'W', 1._wp )
195      CALL lbc_lnk( avm , 'W', 1._wp )
196      CALL lbc_lnk( avmu, 'U', 1._wp ) 
197      CALL lbc_lnk( avmv, 'V', 1._wp )
198
199      IF(ln_ctl) THEN
200         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
201         CALL prt_ctl(tab3d_1=avmu, clinfo1=' ddm  - u: ', mask1=umask, &
202            &         tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk)
203      ENDIF
204      !
205      IF( wrk_not_released(2, 1,2,3,4,5) )   CALL ctl_stop('zdf_ddm: Release of workspace arrays failed')
206      !
207   END SUBROUTINE zdf_ddm
208   
209   
210   SUBROUTINE zdf_ddm_init
211      !!----------------------------------------------------------------------
212      !!                  ***  ROUTINE zdf_ddm_init  ***
213      !!
214      !! ** Purpose :   Initialization of double diffusion mixing scheme
215      !!
216      !! ** Method  :   Read the namzdf_ddm namelist and check the parameter values
217      !!              called by zdf_ddm at the first timestep (nit000)
218      !!----------------------------------------------------------------------
219      NAMELIST/namzdf_ddm/ rn_avts, rn_hsbfr
220      !!----------------------------------------------------------------------
221      !
222      REWIND( numnam )                ! Read Namelist namzdf_ddm : double diffusion mixing scheme
223      READ  ( numnam, namzdf_ddm )
224      !
225      IF(lwp) THEN                    ! Parameter print
226         WRITE(numout,*)
227         WRITE(numout,*) 'zdf_ddm : double diffusive mixing'
228         WRITE(numout,*) '~~~~~~~'
229         WRITE(numout,*) '   Namelist namzdf_ddm : set dd mixing parameter'
230         WRITE(numout,*) '      maximum avs for dd mixing      rn_avts   = ', rn_avts
231         WRITE(numout,*) '      heat/salt buoyancy flux ratio  rn_hsbfr  = ', rn_hsbfr
232      ENDIF
233      !
234      !                              ! allocate zdfddm arrays
235      IF( zdf_ddm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ddm_init : unable to allocate arrays' )
236      !
237   END SUBROUTINE zdf_ddm_init
238
239#else
240   !!----------------------------------------------------------------------
241   !!   Default option :          Dummy module          No double diffusion
242   !!----------------------------------------------------------------------
243   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfddm = .FALSE.   !: double diffusion flag
244CONTAINS
245   SUBROUTINE zdf_ddm( kt )           ! Dummy routine
246      WRITE(*,*) 'zdf_ddm: You should not have seen this print! error?', kt
247   END SUBROUTINE zdf_ddm
248   SUBROUTINE zdf_ddm_init            ! Dummy routine
249   END SUBROUTINE zdf_ddm_init
250#endif
251
252   !!======================================================================
253END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.