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.
ldftra_smag.F90 in branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_smag.F90 @ 4616

Last change on this file since 4616 was 4616, checked in by gm, 10 years ago

#1260 : see the associated wiki page for explanation
File size: 8.1 KB
Line 
1MODULE ldftra_smag
2   !!======================================================================
3   !!                     ***  MODULE  ldftrasmag  ***
4   !! Ocean physics:  variable eddy induced velocity coefficients
5   !!======================================================================
6   !!  last modified : Maria Luneva, October 2012
7   !!----------------------------------------------------------------------
8#if   defined key_traldf_smag
9   !!----------------------------------------------------------------------
10   !!   'key_traldf_smag'      and           smagorinsky  diffusivity
11   !!----------------------------------------------------------------------
12   !!   ldf_tra_smag  : compute the smagorinski eddy coefficients
13   !!----------------------------------------------------------------------
14   USE oce            ! ocean dynamics and tracers
15   USE dom_oce        ! ocean space and time domain
16   USE sbc_oce        ! surface boundary condition: ocean
17   USE sbcrnf         ! river runoffs
18   USE ldftra         ! ocean tracer   lateral physics
19   USE phycst         ! physical constants
20   USE ldfslp         ! iso-neutral slopes
21   !
22   USE in_out_manager ! I/O manager
23   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
24   USE prtctl         ! Print control
25   USE iom            !
26   USE ioipsl         !
27   USE wrk_nemo       !
28   
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC ldf_tra_smag               ! routine called by step.F90
33   
34   !! * Substitutions
35#  include "domzgr_substitute.h90"
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 3.6 , LOCEAN-IPSL (2014)
39   !! $Id: ldf_tra_smag.F90 1482 2010-06-13 15:28:06Z  $
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE ldf_tra_smag( kt )
45      !!----------------------------------------------------------------------
46      !!                  ***  ROUTINE ldf_tra_smag  ***
47      !!                   
48      !! ** Purpose :   initializations of the horizontal ocean physics
49      !!
50      !! ** Method  :   3D eddy viscosity coef.
51      !!    M.Griffies, R.Hallberg AMS, 2000
52      !!  for laplacian:
53      !!   Asmag=(C/pi)^2*dx*dy sqrt(D^2), C=1 for tracers, C=3-4 for viscosity
54      !!   D^2= rm_smsh*(du/dx-dv/dy)^2+(dv/dx+du/dy)^2 for Cartesian coordinates
55      !!  IF rm_smsh = 0 , only shear is used, recommended for tidal flows
56      !!  in general case du/dx ==> e2 d(u/e2)/dx;  du/dy ==> e1 d(u/e1)/dy;
57      !!                  dv/dx ==> e2 d(v/e2)/dx;  dv/dy ==> e1 d(v/e1)/dy
58      !!  for bilaplacian: now this option is deleted as unstable or non-conservative
59      !!   - delta{Bsmag (delta(T)} = -Bsmag* delta{delta(T)} - delta(Bsmag)*delta( T )
60      !!  second term is of arbitrary sign on the edge of fronts and can induce instability
61      !!   Bsmag=Asmag*dx*dy/8
62      !! 
63      !!       laplacian operator   : ahm1, ahm2 defined at T- and F-points
64      !!                              ahm3, ahm4 never used
65      !!       bilaplacian operator : ahm1, ahm2 never used
66      !!                           :  ahm3, ahm4 defined at U- and V-points
67      !!       ??? explanation of the default is missing
68      !!----------------------------------------------------------------------
69      INTEGER, INTENT( in ) ::   kt   ! ocean time-step inedx
70      !
71      INTEGER  ::   ji, jj, jk               ! dummy loop indices
72      REAL(wp) ::   zdeltau, zhsmag          ! local scalars
73      REAL(wp) ::   zdeltav, zsmsh , zcoef   !   -      -
74      REAL(wp), POINTER , DIMENSION (:,:) ::   zux, zvx , zuy , zvy 
75      REAL(wp), POINTER , DIMENSION (:,:) ::   zue1, zue2 , zve1 , zve2 
76     
77      CALL wrk_alloc (jpi,jpj,zux, zvx , zuy , zvy )
78      CALL wrk_alloc (jpi,jpj,zue1, zue2 , zve1 , zve2 )
79      !!----------------------------------------------------------------------
80      !
81      IF(  kt == nit000 ) THEN
82         IF(lwp) WRITE(numout,*)
83         IF(lwp) WRITE(numout,*) ' ldf_tra_smag : 3D eddy smagorinsky diffusivity '
84         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~   --  '
85      ENDIF
86
87      zhsmag = rn_chsmag
88      zsmsh  = rn_smsh 
89      zux(:,:) = 0._wp   ;   zuy(:,:) = 0._wp   ;   zvx(:,:) = 0._wp   ;   zvy(:,:) = 0._wp
90
91      ! -------------------
92      ahtt(:,:,:) = rn_aht_0
93      IF( ln_traldf_bilap ) THEN
94         IF( lwp .AND. kt == nit000) WRITE(numout,* ) 'ldf_tra_smag :no bilaplacian Smagorinsky diffusivity'
95         IF( lwp .AND. kt == nit000) WRITE(numout,* ) 'ldf_tra_smag :bilaplacian diffusivity set to constant' 
96      ENDIF
97
98
99
100      ! harmonic operator   (U-, V-, W-points)
101      ! -----------------
102      ahtu(:,:,:) = rn_aht_0                  ! set ahtu , ahtv at u- and v-points,
103      ahtv(:,:,:) = rn_aht_0                  ! and ahtw at w-point
104     
105      IF( ln_traldf_lap ) THEN
106         !
107         DO jk = 1 , jpkm1
108            zue2(:,:) = un(:,:,jk) * r1_e2u(:,:)          !!gm  for stability reason use of before instead of now here !!!!
109            zve1(:,:) = vn(:,:,jk) * r1_e1v(:,:)
110            zue1(:,:) = un(:,:,jk) * r1_e1u(:,:)
111            zve2(:,:) = vn(:,:,jk) * r1_e2v(:,:)
112            !
113            DO jj = 2, jpj                               !!gm  multiplication by tmask useless as un, vn maked field !
114               DO ji= 2, jpi
115                  zux(ji,jj) = ( zue2(ji,jj) - zue2(ji-1,jj  ) ) * r1_e1e2t(ji,jj) * tmask(ji,jj,jk) * zsmsh 
116                  zvy(ji,jj) = ( zve1(ji,jj) - zve1(ji  ,jj-1) ) * r1_e1e2t(ji,jj) * tmask(ji,jj,jk) * zsmsh 
117               END DO
118            END DO
119            !
120            DO jj = 1, jpjm1
121               DO ji = 1, jpim1
122               zuy(ji,jj) = ( zue1(ji  ,jj+1) - zue1(ji,jj) ) * r1_e2f(ji,jj) *e1f(ji,jj) * fmask(ji,jj,jk)
123               zvx(ji,jj) = ( zve2(ji+1,jj  ) - zve2(ji,jj) ) * r1_e1f(ji,jj) *e2f(ji,jj) * fmask(ji,jj,jk)
124               END DO
125            END DO
126            !
127            DO jj = 2, jpjm1
128               DO ji = 2, jpim1
129                  zdeltau = 2._wp / ( e1u(ji,jj)**(-2) + e2u(ji,jj)**(-2) )
130                  zdeltav = 2._wp / ( e1v(ji,jj)**(-2) + e2v(ji,jj)**(-2) )
131                  !
132                  ahtu(ji,jj,jk) = MAX( rn_aht_0 , (zhsmag/rpi)**2*zdeltau*                                    &
133                     &            SQRT(  0.25_wp*( zux(ji,jj)+zux(ji+1,jj  )-zvy(ji,jj)-zvy(ji+1,jj  ) )**2    &
134                     &                 + 0.25_wp*( zuy(ji,jj)+zuy(ji  ,jj-1)+zvx(ji,jj)+zvx(ji  ,jj-1) )**2  )   )
135                     !
136                  ahtv(ji,jj,jk) = MAX( rn_aht_0 ,  (zhsmag/rpi)**2*zdeltav*                                   &
137                     &             SQRT(  0.25_wp*( zux(ji,jj)+zux(ji  ,jj+1)-zvy(ji  ,jj)-zvy(ji,jj+1) )**2     &
138                     &                  + 0.25_wp*( zuy(ji,jj)+zuy(ji-1,jj  )+zvx(ji-1,jj)+zvx(ji,jj  )  )**2  )   )
139                     !
140                  ! stability criteria: aht<delta**2/(4*dt)   dt=2*rdt , positiveness require aht<delta**2/(8*dt)
141                  ahtu(ji,jj,jk) = MIN( ahtu(ji,jj,jk) , zdeltau / (16*rdt) , rn_aht_m )
142                  ahtv(ji,jj,jk) = MIN( ahtv(ji,jj,jk) , zdeltav / (16*rdt) , rn_aht_m )
143               END DO
144            END DO
145         END DO
146      ENDIF
147      ahtu(:,:,jpk) = ahtu(:,:,jpkm1)
148      ahtv(:,:,jpk) = ahtv(:,:,jpkm1)
149      CALL lbc_lnk( ahtu, 'U', 1. )   ! Lateral boundary conditions
150      CALL lbc_lnk( ahtv, 'V', 1. )
151      !
152      CALL wrk_dealloc ( jpi,jpj,zux, zvx , zuy , zvy     )
153      CALL wrk_dealloc ( jpi,jpj,zue1, zue2 , zve1 , zve2 )
154      !
155   END SUBROUTINE ldf_tra_smag
156
157#else
158   !!----------------------------------------------------------------------
159   !!   Default option                                         Dummy module
160   !!----------------------------------------------------------------------
161   CONTAINS
162   SUBROUTINE ldf_tra_smag( kt )       ! Empty routine
163      WRITE(*,*) 'ldf_dyn_smag: You should not have seen this print! error? check keys ldf:c3d+smag', kt
164   END SUBROUTINE ldf_tra_smag
165#endif
166
167   !!======================================================================
168END MODULE ldftra_smag
Note: See TracBrowser for help on using the repository browser.