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

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

Last change on this file since 5021 was 4990, checked in by timgraham, 9 years ago

Merged branches/2014/dev_MERGE_2014 back onto the trunk as follows:

In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
1 conflict in LIM_SRC_3/limdiahsb.F90
Resolved by keeping the version from dev_MERGE_2014 branch
and commited at r4989

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2014/dev_MERGE_2014
to merge the branch into the trunk - no conflicts at this stage.

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