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 @ 7424

Last change on this file since 7424 was 2931, checked in by hliu, 13 years ago
File size: 10.7 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      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
81      USE wrk_nemo, ONLY:   zux => wrk_2d_1  , zuy => wrk_2d_2  , zvx => wrk_2d_3  , zvy => wrk_2d_4   ! 2D workspace
82      USE wrk_nemo, ONLY:   zue1 => wrk_2d_5 , zue2 => wrk_2d_6 , zve1 => wrk_2d_7 , zve2 => wrk_2d_8   ! 2D workspac
83      INTEGER, INTENT( in ) ::   kt                             ! ocean time-step inedx
84      !! * Arguments
85      INTEGER               :: ji,jj,jk
86
87      REAL (wp)             ::    rdeltaxy, rdeltau, rdeltav    ! temporary scalars
88
89      !!----------------------------------------------------------------------
90      IF(  kt == nit000 ) THEN
91      IF( lk_traldf_eiv ) THEN
92         IF(lwp) WRITE(numout,*)
93         IF(lwp) WRITE(numout,*) ' ldf_tra_smag : 3D eddy smagorinsky diffusivity and eddy induced velocity coefficients'
94         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~   --  '
95         IF(lwp) WRITE(numout,*) '               Coefficients are computed'
96         IF(lwp) WRITE(numout,*)
97      ELSE
98         IF(lwp) WRITE(numout,*)
99         IF(lwp) WRITE(numout,*) ' ldf_tra_smag : 3D eddy smagorinsky diffusivity coefficient'
100         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~   -- '
101         IF(lwp) WRITE(numout,*) '               Coefficients are computed'
102         IF(lwp) WRITE(numout,*)
103      ENDIF
104      ENDIF
105
106      IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN
107         CALL ctl_stop('tra_ldf_smag : requested workspace arrays unavailable')   ;   RETURN
108      ENDIF
109     
110      zux(:,:)=0._wp ; zuy(:,:)=0._wp ; zvx(:,:)=0._wp ; zvy(:,:)=0._wp
111
112      ! biharmonic operator   (T-point)
113      ! -------------------
114     
115      ahtt(:,:,:) = rn_aht_0              ! set ahtt at T-point (here no space variation)
116        IF( ln_traldf_bilap ) THEN
117         ! define ahm1 and ahm2 at the right grid point position
118         ! (USER: modify ahm1 and ahm2 following your desiderata)
119         DO jk=1,jpk
120           zue2(:,:)=un(:,:,jk)/e2u(:,:)
121           zve1(:,:)=vn(:,:,jk)/e1v(:,:)
122           zue1(:,:)=un(:,:,jk)/e1u(:,:)
123           zve2(:,:)=vn(:,:,jk)/e2v(:,:)
124
125
126           DO jj=2,jpj
127            DO ji=2,jpi
128            zux(ji,jj)=(zue2(ji,jj)-zue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
129            zvy(ji,jj)=(zve1(ji,jj)-zve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk)
130            ENDDO
131           ENDDO
132
133           DO jj=1,jpjm1
134            DO ji=1,jpim1
135            zuy(ji,jj)=(zue1(ji,jj+1)-zue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
136            zvx(ji,jj)=(zve2(ji+1,jj)-zve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
137            ENDDO
138           ENDDO
139
140
141          DO jj=2,jpjm1
142           DO ji=2,jpim1
143            rdeltaxy=2./(e1t(ji,jj)**(-2)+e2t(ji,jj)**(-2))
144             ahtt(ji,jj,jk)=(rn_chsmag/3.14)**2/8*rdeltaxy**2*                             &
145                            sqrt( (zux(ji,jj)-zvy(ji,jj))**2+                              &
146                     0.0625*(zuy(ji,jj)+zuy(ji,jj-1)+zuy(ji-1,jj)+zuy(ji-1,jj-1)+          &
147                        zvx(ji,jj)+zvx(ji,jj-1)+zvx(ji-1,jj)+zvx(ji-1,jj-1) )**2)
148   
149              !!! stabulity criteria: abs(aht)<delta**4/(4*8*dt)   dt=2*rdt
150             ahtt(ji,jj,jk)=MIN( ahtt(ji,jj,jk) , rn_aht_0 )     
151             ahtt(ji,jj,jk)=MAX( ahtt(ji,jj,jk) ,- rdeltaxy**4/(8*16*rdt) ) 
152            ENDDO
153           ENDDO
154         ENDDO
155            ahtt(:,:,jpk) = ahtt(:,:,jpkm1)
156
157      CALL lbc_lnk( ahtt, 'T', 1. )   ! Lateral boundary conditions on ( ahtt )
158IF(  kt == nit000 ) THEN
159
160      IF(lwp ) THEN                    ! Control print
161         WRITE(numout,*)
162         WRITE(numout,*) 'inildf: ahtt at k = 1'
163         CALL prihre( ahtt(:,:,1), jpi, jpj, 1, jpi, 1,   &
164            &                                1, jpj, 1, 1.e-1, numout )
165
166      ENDIF
167
168      ENDIF
169ENDIF
170
171
172
173      ! harmonic operator   (U-, V-, W-points)
174      ! -----------------
175
176      ahtu(:,:,:) = rn_aht_0              ! set ahtu = ahtv at u- and v-points,
177      ahtv(:,:,:) = rn_aht_0              ! and ahtw at w-point
178      ahtw(:,:,:) = rn_aht_0              ! (here example: no space variation)
179
180      IF( ln_traldf_lap ) THEN
181         DO jk=1,jpk
182           zue2(:,:)=un(:,:,jk)/e2u(:,:)
183           zve1(:,:)=vn(:,:,jk)/e1v(:,:)
184           zue1(:,:)=un(:,:,jk)/e1u(:,:)
185           zve2(:,:)=vn(:,:,jk)/e2v(:,:)
186
187
188           DO jj=2,jpj
189            DO ji=2,jpi
190            zux(ji,jj)=(zue2(ji,jj)-zue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
191            zvy(ji,jj)=(zve1(ji,jj)-zve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk)
192            ENDDO
193           ENDDO
194
195           DO jj=1,jpjm1
196            DO ji=1,jpim1
197            zuy(ji,jj)=(zue1(ji,jj+1)-zue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
198            zvx(ji,jj)=(zve2(ji+1,jj)-zve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
199            ENDDO
200           ENDDO
201
202
203          DO jj=2,jpjm1
204           DO ji=2,jpim1
205           rdeltau=2./(e1u(ji,jj)**(-2)+e2u(ji,jj)**(-2))
206           rdeltav=2./(e1v(ji,jj)**(-2)+e2v(ji,jj)**(-2))
207
208           ahtu(ji,jj,jk)=MAX(rn_aht_0 , (rn_chsmag/3.14)**2*rdeltau*                       &
209                         sqrt(0.25*(zux(ji,jj)+zux(ji+1,jj)-zvy(ji,jj)-zvy(ji+1,jj))**2+    &
210                              0.25*(zuy(ji,jj)+zuy(ji,jj-1)+zvx(ji,jj)+zvx(ji,jj-1))**2))
211
212           ahtv(ji,jj,jk)=MAX(rn_aht_0 , (rn_chsmag/3.14)**2*rdeltav*                       &
213                         sqrt(0.25*(zux(ji,jj)+zux(ji,jj+1)-zvy(ji,jj)-zvy(ji,jj+1))**2+    &
214                              0.25*(zuy(ji,jj)+zuy(ji-1,jj)+zvx(ji-1,jj)+zvx(ji,jj))**2))
215
216
217         !!! stabulity criteria: aht<delta**2/(4*dt)   dt=2*rdt
218
219             ahtu(ji,jj,jk)=MIN(ahtu(ji,jj,jk),rdeltau**2/(16*rdt) )
220             ahtv(ji,jj,jk)=MIN(ahtv(ji,jj,jk),rdeltav**2/(16*rdt) )
221
222            ENDDO
223           ENDDO
224         ENDDO
225        ENDIF
226            ahtu(:,:,jpk) = ahtu(:,:,jpkm1)
227            ahtv(:,:,jpk) = ahtv(:,:,jpkm1)
228        CALL lbc_lnk( ahtu, 'U', 1. )   ! Lateral boundary conditions
229        CALL lbc_lnk( ahtv, 'V', 1. )
230
231IF(  kt == nit000 ) THEN
232
233      IF(lwp ) THEN                    ! Control print
234         WRITE(numout,*)
235         WRITE(numout,*) 'inildf: ahtu at k = 1'
236         CALL prihre( ahtu(:,:,1), jpi, jpj, 1, jpi, 1,   &
237            &                                1, jpj, 1, 1.e-1, numout )
238         WRITE(numout,*)
239         WRITE(numout,*) 'inildf: ahtv at k = 1'
240         CALL prihre( ahtv(:,:,1), jpi, jpj, 1, jpi, 1,   &
241            &                                1, jpj, 1, 1.e-1, numout )
242         WRITE(numout,*)
243         WRITE(numout,*) 'inildf: ahtw at k = 1'
244         CALL prihre( ahtw(:,:,1), jpi, jpj, 1, jpi, 1,   &
245            &                                1, jpj, 1, 1.e-1, numout )
246      ENDIF
247ENDIF
248
249       IF( wrk_not_released(2, 1,2,3,4,5,6,7,8)) CALL ctl_stop('tra_ldf_smag: failed to release workspace arrays')
250
251END SUBROUTINE ldf_tra_smag
252#else
253   !!----------------------------------------------------------------------
254   !!   Default option                                         Dummy module
255   !!----------------------------------------------------------------------
256CONTAINS
257   SUBROUTINE ldf_tra_smag( kt )       ! Empty routine
258      WRITE(*,*) 'ldf_dyn_smag: You should not have seen this print! error? check keys ldf:c3d+smag', kt
259   END SUBROUTINE ldf_tra_smag
260#endif
261
262END MODULE ldftra_smag
Note: See TracBrowser for help on using the repository browser.