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

source: branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_smag.F90 @ 2887

Last change on this file since 2887 was 2887, checked in by hliu, 13 years ago

addition and modification of files for Smagorinsky method. for Maria Luneva
File size: 10.5 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
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: ldf_tra_smag.F90 1482 2010-06-13 15:28:06Z  $
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      !! for bilaplacian:
64      !!   Bsmag=Asmag*dx*dy/8
65      !!   D^2=(du/dx-dv/dy)^2+(dv/dx+du/dy)^2 for Cartesian coordinates
66      !!  in general case du/dx ==> e2 d(u/e2)/dx;  du/dy ==> e1 d(u/e1)/dy;
67      !!                  dv/dx ==> e2 d(v/e2)/dx;  dv/dy ==> e1 d(v/e1)/dy
68      !! 
69      !!       laplacian operator   : ahm1, ahm2 defined at T- and F-points
70      !!                              ahm3, ahm4 never used
71      !!       bilaplacian operator : ahm1, ahm2 never used
72      !!                           :  ahm3, ahm4 defined at U- and V-points
73      !!       ??? explanation of the default is missing
74      !!  last modified : Maria Luneva, September 2011
75      !!----------------------------------------------------------------------
76      !!
77      !!----------------------------------------------------------------------
78      !! * Modules used
79      USE ioipsl
80      INTEGER, INTENT( in ) ::   kt     ! ocean time-step inedx
81      !! * Arguments
82      INTEGER               :: ji,jj,jk
83      REAL (wp), DIMENSION (:,:), ALLOCATABLE::    ux,uy,vx,vy     ! local velocity derivatives
84      REAL (wp), DIMENSION (:,:), ALLOCATABLE::    ue1,ue2,ve1,ve2 ! local variables
85      REAL (wp)                     ::    deltaxy, deltau, deltav  ! local step
86
87      !!----------------------------------------------------------------------
88      IF(  kt == nit000 ) THEN
89      IF( lk_traldf_eiv ) THEN
90         IF(lwp) WRITE(numout,*)
91         IF(lwp) WRITE(numout,*) ' ldf_tra_smag : 3D eddy smagorinsky diffusivity and eddy induced velocity coefficients'
92         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~   --  '
93         IF(lwp) WRITE(numout,*) '               Coefficients are computed'
94         IF(lwp) WRITE(numout,*)
95      ELSE
96         IF(lwp) WRITE(numout,*)
97         IF(lwp) WRITE(numout,*) ' ldf_tra_smag : 3D eddy smagorinsky diffusivity coefficient'
98         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~   -- '
99         IF(lwp) WRITE(numout,*) '               Coefficients are computed'
100         IF(lwp) WRITE(numout,*)
101      ENDIF
102      ENDIF
103
104      ALLOCATE(ux(jpi,jpj)) ; ux(:,:)=0._wp
105      ALLOCATE(uy(jpi,jpj)) ; uy(:,:)=0._wp
106      ALLOCATE(vx(jpi,jpj)) ; vx(:,:)=0._wp
107      ALLOCATE(vy(jpi,jpj)) ; vy(:,:)=0._wp
108      ALLOCATE(ue1(jpi,jpj)); ALLOCATE(ue2(jpi,jpj)) 
109      ALLOCATE(ve1(jpi,jpj)); ALLOCATE(ve2(jpi,jpj)) 
110
111      ! biharmonic operator   (T-point)
112      ! -------------------
113     
114      ahtt(:,:,:) = aht0              ! set ahtt at T-point (here no space variation)
115        IF( ln_traldf_bilap ) THEN
116         ! define ahm1 and ahm2 at the right grid point position
117         ! (USER: modify ahm1 and ahm2 following your desiderata)
118         DO jk=1,jpk
119           ue2(:,:)=un(:,:,jk)/e2u(:,:)
120           ve1(:,:)=vn(:,:,jk)/e1v(:,:)
121           ue1(:,:)=un(:,:,jk)/e1u(:,:)
122           ve2(:,:)=vn(:,:,jk)/e2v(:,:)
123
124
125           DO jj=2,jpj
126            DO ji=2,jpi
127            ux(ji,jj)=(ue2(ji,jj)-ue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
128            vy(ji,jj)=(ve1(ji,jj)-ve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk)
129            ENDDO
130           ENDDO
131
132           DO jj=1,jpjm1
133            DO ji=1,jpim1
134            uy(ji,jj)=(ue1(ji,jj+1)-ue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
135            vx(ji,jj)=(ve2(ji+1,jj)-ve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
136            ENDDO
137           ENDDO
138
139
140          DO jj=2,jpjm1
141           DO ji=2,jpim1
142            deltaxy=2./(e1t(ji,jj)**(-2)+e2t(ji,jj)**(-2))
143             ahtt(ji,jj,jk)=(chsmag/3.14)**2/8*deltaxy**2*                      &
144                            sqrt( (ux(ji,jj)-vy(ji,jj))**2+                              &
145                     0.0625*(uy(ji,jj)+uy(ji,jj-1)+uy(ji-1,jj)+uy(ji-1,jj-1)+            &
146                        vx(ji,jj)+vx(ji,jj-1)+vx(ji-1,jj)+vx(ji-1,jj-1) )**2)
147   
148              !!! stabulity criteria: abs(aht)<delta**4/(4*8*dt)   dt=2*rdt
149             ahtt(ji,jj,jk)=MIN( ahtt(ji,jj,jk) , aht0 )     
150             ahtt(ji,jj,jk)=MAX( ahtt(ji,jj,jk) ,- deltaxy**4/(8*16*rdt) ) 
151            ENDDO
152           ENDDO
153         ENDDO
154            ahtt(:,:,jpk) = ahtt(:,:,jpkm1)
155
156      CALL lbc_lnk( ahtt, 'T', 1. )   ! Lateral boundary conditions on ( ahtt )
157IF(  kt == nit000 ) THEN
158
159      IF(lwp ) THEN                    ! Control print
160         WRITE(numout,*)
161         WRITE(numout,*) 'inildf: ahtt at k = 1'
162         CALL prihre( ahtt(:,:,1), jpi, jpj, 1, jpi, 1,   &
163            &                                1, jpj, 1, 1.e-1, numout )
164
165      ENDIF
166
167      ENDIF
168ENDIF
169
170
171
172      ! harmonic operator   (U-, V-, W-points)
173      ! -----------------
174
175      ahtu(:,:,:) = aht0              ! set ahtu = ahtv at u- and v-points,
176      ahtv(:,:,:) = aht0              ! and ahtw at w-point
177      ahtw(:,:,:) = aht0              ! (here example: no space variation)
178
179      IF( ln_traldf_lap ) THEN
180         DO jk=1,jpk
181           ue2(:,:)=un(:,:,jk)/e2u(:,:)
182           ve1(:,:)=vn(:,:,jk)/e1v(:,:)
183           ue1(:,:)=un(:,:,jk)/e1u(:,:)
184           ve2(:,:)=vn(:,:,jk)/e2v(:,:)
185
186
187           DO jj=2,jpj
188            DO ji=2,jpi
189            ux(ji,jj)=(ue2(ji,jj)-ue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
190            vy(ji,jj)=(ve1(ji,jj)-ve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk)
191            ENDDO
192           ENDDO
193
194           DO jj=1,jpjm1
195            DO ji=1,jpim1
196            uy(ji,jj)=(ue1(ji,jj+1)-ue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
197            vx(ji,jj)=(ve2(ji+1,jj)-ve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
198            ENDDO
199           ENDDO
200
201
202          DO jj=2,jpjm1
203           DO ji=2,jpim1
204           deltau=2./(e1u(ji,jj)**(-2)+e2u(ji,jj)**(-2))
205           deltav=2./(e1v(ji,jj)**(-2)+e2v(ji,jj)**(-2))
206
207           ahtu(ji,jj,jk)=MAX(aht0 , (chsmag/3.14)**2*deltau*            &
208                         sqrt(0.25*(ux(ji,jj)+ux(ji+1,jj)-vy(ji,jj)-vy(ji+1,jj))**2+    &
209                              0.25*(uy(ji,jj)+uy(ji,jj-1)+vx(ji,jj)+vx(ji,jj-1))**2))
210
211           ahtv(ji,jj,jk)=MAX(aht0 , (chsmag/3.14)**2*deltav*            &
212                         sqrt(0.25*(ux(ji,jj)+ux(ji,jj+1)-vy(ji,jj)-vy(ji,jj+1))**2+    &
213                              0.25*(uy(ji,jj)+uy(ji-1,jj)+vx(ji-1,jj)+vx(ji,jj))**2))
214
215
216         !!! stabulity criteria: aht<delta**2/(4*dt)   dt=2*rdt
217
218             ahtu(ji,jj,jk)=MIN(ahtu(ji,jj,jk),deltau**2/(16*rdt) )
219             ahtv(ji,jj,jk)=MIN(ahtv(ji,jj,jk),deltav**2/(16*rdt) )
220
221            ENDDO
222           ENDDO
223         ENDDO
224        ENDIF
225            ahtu(:,:,jpk) = ahtu(:,:,jpkm1)
226            ahtv(:,:,jpk) = ahtv(:,:,jpkm1)
227        CALL lbc_lnk( ahtu, 'U', 1. )   ! Lateral boundary conditions
228        CALL lbc_lnk( ahtv, 'V', 1. )
229
230
231          DEALLOCATE( ux ) ;DEALLOCATE( uy );DEALLOCATE( vx );DEALLOCATE( vy )     
232          DEALLOCATE( ue1 ) ;DEALLOCATE( ue2 );DEALLOCATE( ve1 );DEALLOCATE( ve2 )
233
234IF(  kt == nit000 ) THEN
235
236      IF(lwp ) THEN                    ! Control print
237         WRITE(numout,*)
238         WRITE(numout,*) 'inildf: ahtu at k = 1'
239         CALL prihre( ahtu(:,:,1), jpi, jpj, 1, jpi, 1,   &
240            &                                1, jpj, 1, 1.e-1, numout )
241         WRITE(numout,*)
242         WRITE(numout,*) 'inildf: ahtv at k = 1'
243         CALL prihre( ahtv(:,:,1), jpi, jpj, 1, jpi, 1,   &
244            &                                1, jpj, 1, 1.e-1, numout )
245         WRITE(numout,*)
246         WRITE(numout,*) 'inildf: ahtw at k = 1'
247         CALL prihre( ahtw(:,:,1), jpi, jpj, 1, jpi, 1,   &
248            &                                1, jpj, 1, 1.e-1, numout )
249      ENDIF
250ENDIF
251
252END SUBROUTINE ldf_tra_smag
253#else
254   !!----------------------------------------------------------------------
255   !!   Default option                                         Dummy module
256   !!----------------------------------------------------------------------
257CONTAINS
258   SUBROUTINE ldf_tra_smag( kt )       ! Empty routine
259      WRITE(*,*) 'ldf_dyn_smag: You should not have seen this print! error? check keys ldf:c3d+smag', kt
260   END SUBROUTINE ldf_tra_smag
261#endif
262
263END MODULE ldftra_smag
Note: See TracBrowser for help on using the repository browser.