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_r4765_CNRS_agrif/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_smag.F90 @ 5581

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

Merged head of trunk into branch

  • Property svn:keywords set to Id
File size: 8.9 KB
Line 
1MODULE ldftra_smag
2   !!======================================================================
3   !!                     ***  MODULE  ldftrasmag  ***
4   !! Ocean physics:  variable eddy induced velocity coefficients
5   !!======================================================================
6#if   defined key_traldf_smag   &&   defined key_traldf_c3d
7   !!----------------------------------------------------------------------
8   !!   'key_traldf_smag'      and           smagorinsky  diffusivity
9   !!   'key_traldf_c3d'                    3D tracer lateral  mixing coef.
10   !!----------------------------------------------------------------------
11   !!   ldf_eiv      : compute the eddy induced velocity coefficients
12   !!----------------------------------------------------------------------
13   !! * Modules used
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_oce      ! ocean tracer   lateral physics
19   USE phycst          ! physical constants
20   USE ldfslp          ! iso-neutral slopes
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 iom
25   USE wrk_nemo
26   IMPLICIT NONE
27   PRIVATE
28
29   !! * Routine accessibility
30   PUBLIC ldf_tra_smag               ! routine called by step.F90
31   !!----------------------------------------------------------------------
32   !!  OPA 9.0 , LOCEAN-IPSL (2005)
33   !! $Id$
34   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
35   !!----------------------------------------------------------------------
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40
41CONTAINS
42
43 
44
45
46
47  !!----------------------------------------------------------------------
48  !!                        ***  ldf_tra_smag.F90  ***
49  !!----------------------------------------------------------------------
50
51
52   SUBROUTINE ldf_tra_smag( kt )
53      !!----------------------------------------------------------------------
54      !!----------------------------------------------------------------------
55      !!                  ***  ROUTINE ldf_tra_smag  ***
56      !!                   
57      !! ** Purpose :   initializations of the horizontal ocean physics
58      !!
59      !! ** Method  :   3D eddy viscosity coef.
60      !!    M.Griffies, R.Hallberg AMS, 2000
61      !!  for laplacian:
62      !!   Asmag=(C/pi)^2*dx*dy sqrt(D^2), C=1 for tracers, C=3-4 for viscosity
63      !!   D^2= rm_smsh*(du/dx-dv/dy)^2+(dv/dx+du/dy)^2 for Cartesian coordinates
64      !!  IF rm_smsh = 0 , only shear is used, recommended for tidal flows
65      !!  in general case du/dx ==> e2 d(u/e2)/dx;  du/dy ==> e1 d(u/e1)/dy;
66      !!                  dv/dx ==> e2 d(v/e2)/dx;  dv/dy ==> e1 d(v/e1)/dy
67      !!  for bilaplacian: now this option is deleted as unstable or non-conservative
68      !!   - delta{Bsmag (delta(T)} = -Bsmag* delta{delta(T)} - delta(Bsmag)*delta( T )
69      !!  second term is of arbitrary sign on the edge of fronts and can induce instability
70      !!   Bsmag=Asmag*dx*dy/8
71      !! 
72      !!       laplacian operator   : ahm1, ahm2 defined at T- and F-points
73      !!                              ahm3, ahm4 never used
74      !!       bilaplacian operator : ahm1, ahm2 never used
75      !!                           :  ahm3, ahm4 defined at U- and V-points
76      !!       ??? explanation of the default is missing
77      !!  last modified : Maria Luneva, October 2012
78      !!----------------------------------------------------------------------
79      !!
80      !!----------------------------------------------------------------------
81      !! * Modules used
82      USE ioipsl
83      REAL ( wp), POINTER , DIMENSION (:,:) :: zux, zvx , zuy , zvy 
84      REAL ( wp), POINTER , DIMENSION (:,:) :: zue1, zue2 , zve1 , zve2 
85      INTEGER, INTENT( in )                 ::   kt                             ! ocean time-step inedx
86      !! * Arguments
87      INTEGER                               :: ji,jj,jk
88
89      REAL (wp)                             ::     zdeltau, zdeltav, zhsmag ,zsmsh    ! temporary scalars
90     
91      CALL wrk_alloc (jpi,jpj,zux, zvx , zuy , zvy )
92      CALL wrk_alloc (jpi,jpj,zue1, zue2 , zve1 , zve2 )
93      !!----------------------------------------------------------------------
94      IF(  kt == nit000 ) THEN
95         IF(lwp) WRITE(numout,*)
96         IF(lwp) WRITE(numout,*) ' ldf_tra_smag : 3D eddy smagorinsky diffusivity '
97         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~   --  '
98         IF(lwp) WRITE(numout,*) '               Coefficients are computed'
99         IF(lwp) WRITE(numout,*)
100         IF(lwp) WRITE(numout,*)
101      ENDIF
102
103      zhsmag = rn_chsmag
104      zsmsh  = rn_smsh 
105      zux(:,:)=0._wp ; zuy(:,:)=0._wp ; zvx(:,:)=0._wp ; zvy(:,:)=0._wp
106
107      ! -------------------
108      ahtt(:,:,:) = rn_aht_0
109       IF( ln_traldf_bilap ) THEN
110        IF( lwp .AND. kt == nit000) WRITE(numout,* )'ldf_tra_smag :no bilaplacian Smagorinsky diffusivity'
111        IF( lwp .AND. kt == nit000) WRITE(numout,* )'ldf_tra_smag :bilaplacian diffusivity set to constant' 
112       ENDIF
113
114
115
116      ! harmonic operator   (U-, V-, W-points)
117      ! -----------------
118
119      ahtu(:,:,:) = rn_aht_0                  ! set ahtu , ahtv at u- and v-points,
120      ahtv(:,:,:) = rn_aht_0                  ! and ahtw at w-point
121      ahtw(:,:,:) = rn_aht_0                  ! (here example: no space variation)
122     
123      IF( ln_traldf_lap ) THEN
124
125         DO jk=1,jpk
126
127           zue2(:,:)=un(:,:,jk)/e2u(:,:)
128           zve1(:,:)=vn(:,:,jk)/e1v(:,:)
129           zue1(:,:)=un(:,:,jk)/e1u(:,:)
130           zve2(:,:)=vn(:,:,jk)/e2v(:,:)
131
132
133           DO jj=2,jpj
134            DO ji=2,jpi
135            zux(ji,jj)=(zue2(ji,jj)-zue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk) * zsmsh 
136            zvy(ji,jj)=(zve1(ji,jj)-zve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk) * zsmsh 
137            ENDDO
138           ENDDO
139
140           DO jj=1,jpjm1
141            DO ji=1,jpim1
142            zuy(ji,jj)=(zue1(ji,jj+1)-zue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
143            zvx(ji,jj)=(zve2(ji+1,jj)-zve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
144            ENDDO
145           ENDDO
146
147
148          DO jj=2,jpjm1
149           DO ji=2,jpim1
150           zdeltau=2._wp/( e1u(ji,jj)**(-2)+e2u(ji,jj)**(-2) )
151           zdeltav=2._wp/( e1v(ji,jj)**(-2)+e2v(ji,jj)**(-2) )
152
153           ahtu(ji,jj,jk)=MAX( rn_aht_0 , (zhsmag/rpi)**2*zdeltau*                                &
154                          SQRT(0.25_wp*( zux(ji,jj)+zux(ji+1,jj)-zvy(ji,jj)-zvy(ji+1,jj) )**2+    &
155                               0.25_wp*( zuy(ji,jj)+zuy(ji,jj-1)+zvx(ji,jj)+zvx(ji,jj-1)  )**2) )
156
157           ahtv(ji,jj,jk)=MAX( rn_aht_0 ,  (zhsmag/rpi)**2*zdeltav*                               &
158                          SQRT(0.25_wp*( zux(ji,jj)+zux(ji,jj+1)-zvy(ji,jj)-zvy(ji,jj+1) )**2+    &
159                               0.25_wp*( zuy(ji,jj)+zuy(ji-1,jj)+zvx(ji-1,jj)+zvx(ji,jj)  )**2) )
160
161
162         !!! stability criteria: aht<delta**2/(4*dt)   dt=2*rdt , positiveness require aht<delta**2/(8*dt)
163             ahtu(ji,jj,jk)=MIN(ahtu(ji,jj,jk),zdeltau/(16*rdt) ,rn_aht_m)
164             ahtv(ji,jj,jk)=MIN(ahtv(ji,jj,jk),zdeltav/(16*rdt) ,rn_aht_m)
165         ! so...
166
167
168            ENDDO
169           ENDDO
170         ENDDO
171        ENDIF
172            ahtu(:,:,jpk) = ahtu(:,:,jpkm1)
173            ahtv(:,:,jpk) = ahtv(:,:,jpkm1)
174        CALL lbc_lnk( ahtu, 'U', 1. )   ! Lateral boundary conditions
175        CALL lbc_lnk( ahtv, 'V', 1. )
176
177IF(  kt == nit000 ) THEN
178
179      IF(lwp ) THEN                    ! Control print
180         WRITE(numout,*)
181         WRITE(numout,*) 'inildf: ahtu at k = 1'
182         CALL prihre( ahtu(:,:,1), jpi, jpj, 1, jpi, 1,   &
183            &                                1, jpj, 1, 1.e-1, numout )
184         WRITE(numout,*)
185         WRITE(numout,*) 'inildf: ahtv at k = 1'
186         CALL prihre( ahtv(:,:,1), jpi, jpj, 1, jpi, 1,   &
187            &                                1, jpj, 1, 1.e-1, numout )
188         WRITE(numout,*)
189         WRITE(numout,*) 'inildf: ahtw at k = 1'
190         CALL prihre( ahtw(:,:,1), jpi, jpj, 1, jpi, 1,   &
191            &                                1, jpj, 1, 1.e-1, numout )
192      ENDIF
193ENDIF
194
195      CALL wrk_dealloc ( jpi,jpj,zux, zvx , zuy , zvy     )
196      CALL wrk_dealloc ( jpi,jpj,zue1, zue2 , zve1 , zve2 )
197
198
199END SUBROUTINE ldf_tra_smag
200#else
201   !!----------------------------------------------------------------------
202   !!   Default option                                         Dummy module
203   !!----------------------------------------------------------------------
204CONTAINS
205   SUBROUTINE ldf_tra_smag( kt )       ! Empty routine
206      WRITE(*,*) 'ldf_dyn_smag: You should not have seen this print! error? check keys ldf:c3d+smag', kt
207   END SUBROUTINE ldf_tra_smag
208#endif
209
210END MODULE ldftra_smag
Note: See TracBrowser for help on using the repository browser.