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

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 9019

Last change on this file since 9019 was 9019, checked in by timgraham, 6 years ago

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