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

source: tags/nemo_dev_x6/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90 @ 8525

Last change on this file since 8525 was 147, checked in by opalod, 20 years ago

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