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 @ 4424

Last change on this file since 4424 was 4424, checked in by trackstand2, 10 years ago

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