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

source: trunk/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90 @ 605

Last change on this file since 605 was 473, checked in by opalod, 18 years ago

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