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

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

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