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 NEMO/branches/2019/fix_sn_cfctl_ticket2328/src/OCE/ZDF – NEMO

source: NEMO/branches/2019/fix_sn_cfctl_ticket2328/src/OCE/ZDF/zdfddm.F90 @ 11872

Last change on this file since 11872 was 11872, checked in by acc, 4 years ago

Branch 2019/fix_sn_cfctl_ticket2328. See #2328. Replacement of ln_ctl and activation of full functionality with
sn_cfctl structure. These changes rename structure components l_mppout and l_mpptop as l_prtctl and l_prttrc
and introduce l_glochk to activate former ln_ctl code in stpctl.F90 to perform global location of min and max
checks. Also added is l_allon which can be used to activate all output (much like the former ln_ctl). If l_allon
is .false. then l_config decides whether or not the suboptions are used.

   sn_cfctl%l_glochk = .FALSE.    ! Range sanity checks are local (F) or global (T). Set T for debugging only
   sn_cfctl%l_allon  = .FALSE.    ! IF T activate all options. If F deactivate all unless l_config is T
   sn_cfctl%l_config = .TRUE.     ! IF .true. then control which reports are written with the remaining options

Note, these changes pass SETTE tests but all references to ln_ctl need to be removed from the sette scripts.
  • Property svn:keywords set to Id
File size: 8.9 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
[9019]10   !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only
[1601]11   !!----------------------------------------------------------------------
[9019]12
[3]13   !!----------------------------------------------------------------------
[9019]14   !!   zdf_ddm       : compute the Kz for salinity
[3]15   !!----------------------------------------------------------------------
[9019]16   USE oce            ! ocean dynamics and tracers variables
17   USE dom_oce        ! ocean space and time domain variables
18   USE zdf_oce        ! ocean vertical physics variables
[4990]19   USE eosbn2         ! equation of state
20   !
[9019]21   USE in_out_manager ! I/O manager
22   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
23   USE prtctl         ! Print control
24   USE lib_mpp        ! MPP library
[3]25
26   IMPLICIT NONE
27   PRIVATE
28
[2528]29   PUBLIC   zdf_ddm       ! called by step.F90
[3]30
31   !! * Substitutions
32#  include "vectopt_loop_substitute.h90"
33   !!----------------------------------------------------------------------
[9598]34   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[2528]35   !! $Id$
[10068]36   !! Software governed by the CeCILL license (see ./LICENSE)
[3]37   !!----------------------------------------------------------------------
38CONTAINS
39
[9019]40   SUBROUTINE zdf_ddm( kt, p_avm, p_avt, p_avs )
[2715]41      !!----------------------------------------------------------------------
[3]42      !!                  ***  ROUTINE zdf_ddm  ***
43      !!                   
44      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
[1601]45      !!              effect of salt fingering and diffusive convection.
[3]46      !!
47      !! ** Method  :   Diapycnal mixing is increased in case of double
48      !!      diffusive mixing (i.e. salt fingering and diffusive layering)
49      !!      following Merryfield et al. (1999). The rate of double diffusive
[4990]50      !!      mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]):
[3]51      !!         * salt fingering (Schmitt 1981):
[4990]52      !!      for R > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (R/rn_hsbfr)^6 )
53      !!      for R > 1 and rn2 > 0 : zavfs = O
54      !!      otherwise                : zavft = 0.7 zavs / R
[3]55      !!         * diffusive layering (Federov 1988):
[4990]56      !!      for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) )
[3]57      !!      otherwise                   : zavdt = 0
[4990]58      !!      for .5 < R < 1 and N^2 > 0 : zavds = zavdt (1.885 R -0.85)
59      !!      for  0 < R <.5 and N^2 > 0 : zavds = zavdt 0.15 R     
[3]60      !!      otherwise                     : zavds = 0
61      !!         * update the eddy diffusivity:
62      !!      avt = avt + zavft + zavdt
63      !!      avs = avs + zavfs + zavds
[9019]64      !!      avm is required to remain at least above avt and avs.
[3]65      !!     
[1601]66      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
[3]67      !!
[1601]68      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
[3]69      !!----------------------------------------------------------------------
[9019]70      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step indexocean time step
71      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm   !  Kz on momentum    (w-points)
72      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avt   !  Kz on temperature (w-points)
73      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avs   !  Kz on salinity    (w-points)
[2715]74      !
[1601]75      INTEGER  ::   ji, jj , jk     ! dummy loop indices
[4990]76      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
77      REAL(wp) ::   zdt, zds
78      REAL(wp) ::   zinr, zrr       !   -      -
79      REAL(wp) ::   zavft, zavfs    !   -      -
80      REAL(wp) ::   zavdt, zavds    !   -      -
[9019]81      REAL(wp), DIMENSION(jpi,jpj) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
[3]82      !!----------------------------------------------------------------------
[3294]83      !
[3]84      !                                                ! ===============
85      DO jk = 2, jpkm1                                 ! Horizontal slab
86         !                                             ! ===============
87         ! Define the mask
88         ! ---------------
[9019]89!!gm  WORK to be done:   change the code from vector optimisation to scalar one.
90!!gm                     ==>>>  test in the loop instead of use of mask arrays
91!!gm                            and many acces in memory
92         
93         DO jj = 1, jpj                !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
[4990]94            DO ji = 1, jpi
[6140]95               zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
[9019]96!!gm please, use e3w_n below
[6140]97                  &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) ) 
[4990]98               !
99               zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  )  &
100                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
101               zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  )  &
102                   &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
103               !
104               zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )
105               zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) 
106               IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
107               zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau
108            END DO
109         END DO
[3]110
[9019]111         DO jj = 1, jpj                !==  indicators  ==!
[3]112            DO ji = 1, jpi
113               ! stability indicator: msks=1 if rn2>0; 0 elsewhere
[2715]114               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
115               ELSE                                       ;   zmsks(ji,jj) = 1._wp
[3]116               ENDIF
[4990]117               ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere           
118               IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp
[2715]119               ELSE                                       ;   zmskf(ji,jj) = 1._wp
[3]120               ENDIF
121               ! diffusive layering indicators:
[4990]122               !     ! mskdl1=1 if 0< R <1; 0 elsewhere
123               IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp
[2715]124               ELSE                                       ;   zmskd1(ji,jj) = 1._wp
[3]125               ENDIF
[4990]126               !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere
127               IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp
[2715]128               ELSE                                       ;   zmskd2(ji,jj) = 1._wp
[3]129               ENDIF
[4990]130               !   mskdl3=1 if 0.5< R <1; 0 elsewhere
131               IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
132               ELSE                                                   ;   zmskd3(ji,jj) = 1._wp
[3]133               ENDIF
134            END DO
135         END DO
136         ! mask zmsk in order to have avt and avs masked
[7753]137         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk)
[3]138
139
140         ! Update avt and avs
141         ! ------------------
142         ! Constant eddy coefficient: reset to the background value
143         DO jj = 1, jpj
144            DO ji = 1, jpi
[4990]145               zinr = 1._wp / zrau(ji,jj)
[3]146               ! salt fingering
[4990]147               zrr = zrau(ji,jj) / rn_hsbfr
[3]148               zrr = zrr * zrr
[1601]149               zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
[1163]150               zavft = 0.7 * zavfs * zinr
[3]151               ! diffusive layering
[1601]152               zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
[4990]153               zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   &
154                  &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
[3]155               ! add to the eddy viscosity coef. previously computed
[9019]156               p_avs(ji,jj,jk) = p_avt(ji,jj,jk) + zavfs + zavds
157               p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zavft + zavdt
158               p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
[3]159            END DO
160         END DO
161         !                                                ! ===============
162      END DO                                              !   End of slab
163      !                                                   ! ===============
[1601]164      !
[11872]165      IF(sn_cfctl%l_prtctl) THEN
[9440]166         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', kdim=jpk)
[49]167      ENDIF
[1601]168      !
[3]169   END SUBROUTINE zdf_ddm
170   
171   !!======================================================================
172END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.