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

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 9366

Last change on this file since 9366 was 9366, checked in by andmirek, 6 years ago

#2050 first version. Compiled OK in moci test suite

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