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

source: branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90 @ 2240

Last change on this file since 2240 was 2240, checked in by cetlod, 13 years ago

Suppression of key_zco everywhere in the code

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