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.
ldfdyn_c3d.h90 in branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90 @ 2436

Last change on this file since 2436 was 2436, checked in by gm, 13 years ago

v3.3beta: suppress obsolete key_orca_lev10
  • Property svn:keywords set to Id
File size: 17.0 KB
Line 
1   !!----------------------------------------------------------------------
2   !!                        ***  ldfdyn_c3d.h90  ***
3   !!----------------------------------------------------------------------
4
5   !!----------------------------------------------------------------------
6   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
7   !! $Id$
8   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!   'key_dynldf_c3d'             3D lateral eddy viscosity coefficients
13   !!----------------------------------------------------------------------
14
15   SUBROUTINE ldf_dyn_c3d( ld_print )
16      !!----------------------------------------------------------------------
17      !!                  ***  ROUTINE ldf_dyn_c3d  ***
18      !!                   
19      !! ** Purpose :   initializations of the horizontal ocean physics
20      !!
21      !! ** Method  :   3D eddy viscosity coef. ( longitude, latitude, depth )
22      !!       laplacian operator   : ahm1, ahm2 defined at T- and F-points
23      !!                              ahm2, ahm4 never used
24      !!       bilaplacian operator : ahm1, ahm2 never used
25      !!                           :  ahm3, ahm4 defined at U- and V-points
26      !!       ??? explanation of the default is missing
27      !!----------------------------------------------------------------------
28      USE ldftra_oce, ONLY : aht0
29      !!
30      LOGICAL, INTENT (in) :: ld_print   ! If true, output arrays on numout
31      !!
32      INTEGER  ::   ji, jj, jk      ! dummy loop indices
33      REAL(wp) ::   &
34         zr = 0.2 ,   &  ! maximum of the reduction factor at the bottom ocean
35         !               ! ( 0 < zr < 1 )
36         zh = 500.,   &  ! depth of at which start the reduction ( > dept(1) )
37         zd_max   ,   &  ! maximum grid spacing over the global domain
38         za00, zc, zd    ! temporary scalars
39      REAL(wp) ::        &
40         zetmax, zefmax, &
41         zeumax, zevmax   
42      REAL(wp), DIMENSION(jpk) ::   zcoef   ! temporary workspace
43      !!----------------------------------------------------------------------
44
45      IF(lwp) WRITE(numout,*)
46      IF(lwp) WRITE(numout,*) 'ldf_dyn_c3d : 3D lateral eddy viscosity coefficient'
47      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
48
49     
50      ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operators
51      ! =================                       whatever its orientation is)
52      IF( ln_dynldf_lap ) THEN
53         ! define ahm1 and ahm2 at the right grid point position
54         ! (USER: modify ahm1 and ahm2 following your desiderata)
55
56         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
57         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
58
59         IF(lwp) WRITE(numout,*) '              laplacian operator: ahm proportional to e1'
60         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
61
62         za00 = ahm0 / zd_max
63
64         IF( ln_dynldf_iso ) THEN
65            IF(lwp) WRITE(numout,*) '              Caution, as implemented now, the isopycnal part of momentum'
66            IF(lwp) WRITE(numout,*) '                 mixing use aht0 as eddy viscosity coefficient. Thus, it is'
67            IF(lwp) WRITE(numout,*) '                 uniform and you must be sure that your ahm is greater than'
68            IF(lwp) WRITE(numout,*) '                 aht0 everywhere in the model domain.'
69         ENDIF
70
71         CALL ldf_zpf( .TRUE. , 1000., 500., 0.25, fsdept(:,:,:), ahm1 )   ! vertical profile
72         CALL ldf_zpf( .TRUE. , 1000., 500., 0.25, fsdept(:,:,:), ahm2 )   ! vertical profile
73         DO jk = 1,jpk
74            DO jj = 1, jpj
75               DO ji = 1, jpi
76                  zetmax = MAX( e1t(ji,jj), e2t(ji,jj) )
77                  zefmax = MAX( e1f(ji,jj), e2f(ji,jj) )
78                  ahm1(ji,jj,jk) = za00 * zetmax * ahm1(ji,jj,jk)
79                  ahm2(ji,jj,jk) = za00 * zefmax * ahm2(ji,jj,jk)
80               END DO
81            END DO
82         END DO
83
84
85         ! Special case for ORCA R2 and R4 configurations (overwrite the value of ahm1 ahm2)
86         ! ==============================================
87         IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) ) THEN
88            IF(lwp) WRITE(numout,*)
89            IF(lwp) WRITE(numout,*) '              ORCA R2 or R4: overwrite the previous definition of ahm'
90            IF(lwp) WRITE(numout,*) '              ============='
91            CALL ldf_dyn_c3d_orca( ld_print )
92         ENDIF
93
94      ENDIF
95     
96      ! Control print
97      IF(lwp .AND. ld_print ) THEN
98         WRITE(numout,*)
99         WRITE(numout,*) '         3D ahm1 array (k=1)'
100         CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
101         WRITE(numout,*)
102         WRITE(numout,*) '         3D ahm2 array (k=1)'
103         CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
104      ENDIF
105
106
107      ! ahm3 and ahm4 at U- and V-points (used for bilaplacian operator
108      ! ================================  whatever its orientation is)
109      ! (USER: modify ahm3 and ahm4 following your desiderata)
110      ! Here: ahm is proportional to the cube of the maximum of the gridspacing
111      !       in the to horizontal direction
112
113      IF( ln_dynldf_bilap ) THEN
114
115         zd_max = MAX( MAXVAL( e1u(:,:) ), MAXVAL( e2u(:,:) ) )
116         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
117
118         IF(lwp) WRITE(numout,*) '              bi-laplacian operator: ahm proportional to e1**3 '
119         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
120
121         za00 = ahm0_blp / ( zd_max * zd_max * zd_max )
122         DO jj = 1, jpj
123            DO ji = 1, jpi
124               zeumax = MAX( e1u(ji,jj), e2u(ji,jj) )
125               zevmax = MAX( e1v(ji,jj), e2v(ji,jj) )
126               ahm3(ji,jj,1) = za00 * zeumax * zeumax * zeumax
127               ahm4(ji,jj,1) = za00 * zevmax * zevmax * zevmax
128            END DO
129         END DO
130
131         zh = MAX( zh, fsdept(1,1,1) )   ! at least the first reach ahm0
132         IF( ln_zco ) THEN               ! z-coordinate, same profile everywhere
133            IF(lwp) WRITE(numout,'(36x," ahm ", 7x)')
134            DO jk = 1, jpk
135               IF( fsdept(1,1,jk) <= zh ) THEN
136                  zcoef(jk) = 1.e0
137               ELSE
138                  zcoef(jk) = 1.e0 + ( zr - 1.e0 )   &
139                     &               * (  1. - EXP( ( fsdept(1,1,jk   ) - zh ) / zh )  )   &
140                     &               / (  1. - EXP( ( fsdept(1,1,jpkm1) - zh ) / zh )  )
141               ENDIF
142               ahm3(:,:,jk) = ahm3(:,:,1) * zcoef(jk)
143               ahm4(:,:,jk) = ahm4(:,:,1) * zcoef(jk)
144               IF(lwp) WRITE(numout,'(34x,E7.2,8x,i3)') zcoef(jk) * ahm0, jk
145            END DO
146         ELSE                            ! partial steps or s-ccordinate
147            zc = MAXVAL( fsdept(:,:,jpkm1) )
148            IF( lk_mpp )   CALL mpp_max( zc )   ! max over the global domain
149
150            zc = 1. / (  1. - EXP( ( zc - zh ) / zh )  )
151            DO jk = 2, jpkm1
152               DO jj = 1, jpj
153                  DO ji = 1, jpi
154                     IF( fsdept(ji,jj,jk) <= zh ) THEN
155                        ahm3(ji,jj,jk) = ahm3(ji,jj,1)
156                        ahm4(ji,jj,jk) = ahm4(ji,jj,1)
157                     ELSE
158                        zd = 1.e0 + ( zr - 1.e0 ) * (  1. - EXP( ( fsdept(ji,jj,jk) - zh ) / zh )  ) * zc
159                        ahm3(ji,jj,jk) = ahm3(ji,jj,1) * zd
160                        ahm4(ji,jj,jk) = ahm4(ji,jj,1) * zd
161                     ENDIF
162                  END DO
163               END DO
164            END DO
165            ahm3(:,:,jpk) = ahm3(:,:,jpkm1)
166            ahm4(:,:,jpk) = ahm4(:,:,jpkm1)
167            IF(lwp) WRITE(numout,'(36x," ahm ", 7x)')
168            DO jk = 1, jpk
169               IF(lwp) WRITE(numout,'(30x,E10.2,8x,i3)') ahm3(1,1,jk), jk
170            END DO
171         ENDIF
172
173         ! Control print
174         IF( lwp .AND. ld_print ) THEN
175            WRITE(numout,*)
176            WRITE(numout,*) 'inildf: ahm3 array at level 1'
177            CALL prihre(ahm3(:,:,1  ),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
178            WRITE(numout,*)
179            WRITE(numout,*) 'inildf: ahm4 array at level 1'
180            CALL prihre(ahm4(:,:,jpk),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
181         ENDIF
182      ENDIF
183
184   END SUBROUTINE ldf_dyn_c3d
185
186
187   SUBROUTINE ldf_dyn_c3d_orca( ld_print )
188      !!----------------------------------------------------------------------
189      !!                  ***  ROUTINE ldf_dyn_c3d  ***
190      !!                   
191      !! ** Purpose :   ORCA R2 an R4 only
192      !!
193      !! ** Method  :   blah blah blah ....
194      !!----------------------------------------------------------------------
195      USE ldftra_oce, ONLY : aht0
196      !!
197      LOGICAL, INTENT (in) :: ld_print   ! If true, output arrays on numout
198      !!
199      INTEGER ::   ji, jj, jk, jn      ! dummy loop indices
200      INTEGER ::   ii0, ii1, ij0, ij1  ! temporary integers
201      INTEGER ::   inum                ! temporary logical unit
202      INTEGER ::   iim, ijm
203      INTEGER ::   ifreq, il1, il2, ij, ii
204      INTEGER, DIMENSION(jpidta, jpjdta) ::   idata
205      INTEGER, DIMENSION(jpi   , jpj   ) ::   icof
206
207      REAL(wp) ::   &
208         zahmeq, zcoff, zcoft, zmsk,   & ! ???
209         zemax, zemin, zeref, zahmm
210      REAL(wp), DIMENSION(jpi,jpj) ::   zahm0
211      REAL(wp), DIMENSION(jpk) ::   zcoef
212
213      CHARACTER (len=15) ::   clexp
214      !!----------------------------------------------------------------------
215
216      IF(lwp) WRITE(numout,*)
217      IF(lwp) WRITE(numout,*) 'ldfdyn_c3d_orca : 3D eddy viscosity coefficient'
218      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
219      IF(lwp) WRITE(numout,*)
220      IF(lwp) WRITE(numout,*) '        orca R2 or R4 ocean model'
221      IF(lwp) WRITE(numout,*) '  reduced in the surface Eq. strip '
222      IF(lwp) WRITE(numout,*)
223
224      ! Read 2d integer array to specify western boundary increase in the
225      ! ===================== equatorial strip (20N-20S) defined at t-points
226
227      CALL ctl_opn( inum, 'ahmcoef', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
228      READ(inum,9101) clexp, iim, ijm
229      READ(inum,'(/)')
230      ifreq = 40
231      il1 = 1
232      DO jn = 1, jpidta/ifreq+1
233         READ(inum,'(/)')
234         il2 = MIN( jpidta, il1+ifreq-1 )
235         READ(inum,9201) ( ii, ji = il1, il2, 5 )
236         READ(inum,'(/)')
237         DO jj = jpjdta, 1, -1
238            READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
239         END DO
240         il1 = il1 + ifreq
241      END DO
242     
243      DO jj = 1, nlcj
244         DO ji = 1, nlci
245            icof(ji,jj) = idata( mig(ji), mjg(jj) )
246         END DO
247      END DO
248      DO jj = nlcj+1, jpj
249         DO ji = 1, nlci
250            icof(ji,jj) = icof(ji,nlcj)
251         END DO
252      END DO
253      DO jj = 1, jpj
254         DO ji = nlci+1, jpi
255            icof(ji,jj) = icof(nlci,jj)
256         END DO
257      END DO
258     
2599101  FORMAT(1x,a15,2i8)
2609201  FORMAT(3x,13(i3,12x))
2619202  FORMAT(i3,41i3)
262     
263      ! Set ahm1 and ahm2
264      ! =================
265     
266      ! define ahm1 and ahm2 at the right grid point position
267      ! (USER: modify ahm1 and ahm2 following your desiderata)
268      ! biharmonic : ahm1 (ahm2) defined at u- (v-) point
269      ! harmonic   : ahm1 (ahm2) defined at t- (f-) point
270     
271      ! first level : as for 2D coefficients
272     
273      ! Decrease ahm to zahmeq m2/s in the tropics
274      ! (from 90 to 20 degre: ahm = constant
275      ! from 20 to  2.5 degre: ahm = decrease in (1-cos)/2
276      ! from  2.5 to  0 degre: ahm = constant
277      ! symmetric in the south hemisphere)
278     
279      IF( jp_cfg == 4 )   THEN
280         zahmeq = 5.0 * aht0
281         zahmm  = min( 160000.0, ahm0)
282         zemax = MAXVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 )
283         zemin = MINVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 )
284         zeref = MAXVAL ( e1t(:,:) * e2t(:,:),   &
285             &   tmask(:,:,1) .GE. 0.5 .AND. ABS(gphit(:,:)) .GT. 50. )
286 
287         DO jj = 1, jpj
288           DO ji = 1, jpi
289              zmsk = e1t(ji,jj) * e2t(ji,jj)
290              IF( abs(gphit(ji,jj)) .LE. 15 ) THEN
291                 zahm0(ji,jj) = ahm0
292              ELSE
293                 IF( zmsk .GE. zeref ) THEN
294                    zahm0(ji,jj) = ahm0
295                 ELSE
296                    zahm0(ji,jj) = zahmm + (ahm0-zahmm)*(1.0 -   &
297                        &          cos((rpi*0.5*(zmsk-zemin)/(zeref-zemin))))
298                 ENDIF
299              ENDIF
300           END DO
301         END DO
302      ENDIF
303
304      IF( jp_cfg == 2 )   THEN
305         zahmeq     = aht0
306         zahmm      = ahm0
307         zahm0(:,:) = ahm0
308      ENDIF
309
310      DO jj = 1, jpj
311         DO ji = 1, jpi
312            IF( ABS(gphif(ji,jj)) >= 20.) THEN
313               ahm2(ji,jj,1) =  zahm0(ji,jj)
314            ELSEIF( ABS(gphif(ji,jj)) <= 2.5) THEN
315               ahm2(ji,jj,1) =  zahmeq
316            ELSE
317               ahm2(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2.   &
318                  &            *(1.-COS( rad*(ABS(gphif(ji,jj))-2.5)*180./17.5 ) )
319            ENDIF
320            IF( ABS(gphit(ji,jj)) >= 20.) THEN
321               ahm1(ji,jj,1) =  zahm0(ji,jj)
322            ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
323               ahm1(ji,jj,1) =  zahmeq
324            ELSE
325               ahm1(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2.   &
326                  &            *(1.-COS( rad*(ABS(gphit(ji,jj))-2.5)*180./17.5 ) )
327            ENDIF
328         END DO
329      END DO
330     
331      ! increase along western boundaries of equatorial strip
332      ! t-point
333      DO jj = 1, jpjm1
334         DO ji = 1, jpim1
335            zcoft = float( icof(ji,jj) ) / 100.
336            ahm1(ji,jj,1) = zcoft * zahm0(ji,jj) + (1.-zcoft) * ahm1(ji,jj,1)
337         END DO
338      END DO
339      ! f-point
340      icof(:,:) = icof(:,:) * tmask(:,:,1)
341      DO jj = 1, jpjm1
342         DO ji = 1, jpim1   ! NO vector opt.
343            zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
344            IF( zmsk == 0. ) THEN
345               zcoff = 1.
346            ELSE
347               zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   &
348                     / (zmsk * 100.)
349            ENDIF
350            ahm2(ji,jj,1) = zcoff * zahm0(ji,jj) + (1.-zcoff) * ahm2(ji,jj,1)
351         END DO
352      END DO
353
354      ! other level: re-increase the coef in the deep ocean
355       DO jk = 1, 21
356         zcoef(jk) = 1._wp
357      END DO
358      zcoef(22) = 2._wp
359      zcoef(23) = 3._wp
360      zcoef(24) = 5._wp
361      zcoef(25) = 7._wp
362      zcoef(26) = 9._wp
363      DO jk = 27, jpk
364         zcoef(jk) = 10._wp
365      END DO
366      !
367      DO jk = 2, jpk
368         ahm1(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm1(:,:,1) )
369         ahm2(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm2(:,:,1) )
370      END DO
371
372      IF( jp_cfg == 4 )   THEN         ! Limit AHM in Gibraltar strait
373         ij0 = 50   ;   ij1 = 53
374         ii0 = 69   ;   ii1 = 71
375         DO jk = 1, jpk
376            ahm1(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) = MIN( zahmm, ahm1(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) )
377            ahm2(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) = MIN( zahmm, ahm2(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) )
378         END DO
379      ENDIF
380      CALL lbc_lnk( ahm1, 'T', 1. )   ! Lateral boundary conditions (unchanged sign)
381      CALL lbc_lnk( ahm2, 'F', 1. )
382
383
384      IF(lwp) THEN                    ! Control print
385         WRITE(numout,*)
386         WRITE(numout,*) '         3D ahm1 array (k=1)'
387         CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
388         WRITE(numout,*)
389         WRITE(numout,*) '         3D ahm2 array (k=1)'
390         CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
391         WRITE(numout,*)
392         WRITE(numout,*) '         3D ahm2 array (k=jpk)'
393         CALL prihre( ahm2(:,:,jpk), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
394      ENDIF
395
396
397      ! Set ahm3 and ahm4
398      ! =================
399
400      ! define ahm3 and ahm4 at the right grid point position
401      ! initialization to a constant value
402      !     (USER: modify ahm3 and ahm4 following your desiderata)
403      !     harmonic isopycnal or geopotential:
404      !                          ahm3 (ahm4) defined at u- (v-) point
405      DO jk = 1, jpk
406         DO jj = 2, jpj
407            DO ji = 2, jpi
408               ahm3(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji  ,jj-1,jk) )
409               ahm4(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji-1,jj  ,jk) )
410            END DO
411         END DO
412      END DO
413      ahm3 ( :, 1, :) = ahm3 ( :, 2, :)
414      ahm4 ( :, 1, :) = ahm4 ( :, 2, :)
415     
416      CALL lbc_lnk( ahm3, 'U', 1. )    ! Lateral boundary conditions (unchanged sign)
417      CALL lbc_lnk( ahm4, 'V', 1. )
418
419      ! Control print
420
421      IF( lwp .AND. ld_print ) THEN
422         WRITE(numout,*)
423         WRITE(numout,*) '         ahm3 array level 1'
424         CALL prihre(ahm3(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
425         WRITE(numout,*)
426         WRITE(numout,*) '         ahm4 array level 1'
427         CALL prihre(ahm4(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
428      ENDIF
429
430   END SUBROUTINE ldf_dyn_c3d_orca
Note: See TracBrowser for help on using the repository browser.