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_c2d.h90 in branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90 @ 4325

Last change on this file since 4325 was 4325, checked in by cbricaud, 10 years ago

dev_MERGE_2013: solve tickets 1192, 1193, 1194
  • Property svn:keywords set to Id
File size: 21.9 KB
Line 
1   !!----------------------------------------------------------------------
2   !!                      ***  ldfdyn_c2d.h90  ***
3   !!----------------------------------------------------------------------
4   !!   ldf_dyn_c2d  : set the lateral viscosity coefficients
5   !!   ldf_dyn_c2d_orca : specific case for orca r2 and r4
6   !!----------------------------------------------------------------------
7
8   !!----------------------------------------------------------------------
9   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
10   !! $Id$
11   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
12   !!----------------------------------------------------------------------
13
14   SUBROUTINE ldf_dyn_c2d( ld_print )
15      !!----------------------------------------------------------------------
16      !!                 ***  ROUTINE ldf_dyn_c2d  ***
17      !!                 
18      !! ** Purpose :   initializations of the horizontal ocean physics
19      !!
20      !! ** Method :
21      !!      2D eddy viscosity coefficients ( longitude, latitude )
22      !!
23      !!       harmonic operator   : ahm1 is defined at t-point
24      !!                             ahm2 is defined at f-point
25      !!           + isopycnal     : ahm3 is defined at u-point
26      !!           or geopotential   ahm4 is defined at v-point
27      !!           iso-model level : ahm3, ahm4 not used
28      !!
29      !!       biharmonic operator : ahm3 is defined at u-point
30      !!                             ahm4 is defined at v-point
31      !!                           : ahm1, ahm2 not used
32      !!
33      !!----------------------------------------------------------------------
34      LOGICAL, INTENT (in) :: ld_print   ! If true, output arrays on numout
35      !
36      INTEGER  ::   ji, jj
37      REAL(wp) ::   za00, zd_max, zetmax, zeumax, zefmax, zevmax
38      !!----------------------------------------------------------------------
39
40      IF(lwp) WRITE(numout,*)
41      IF(lwp) WRITE(numout,*) 'ldf_dyn_c2d : 2d lateral eddy viscosity coefficient'
42      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
43
44      ! harmonic operator (ahm1, 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         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
51         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
52
53         IF(lwp) WRITE(numout,*) '              laplacian operator: ahm proportional to e1'
54         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
55
56         za00 = ahm0 / zd_max
57         DO jj = 1, jpj
58            DO ji = 1, jpi
59               zetmax = MAX( e1t(ji,jj), e2t(ji,jj) )
60               zefmax = MAX( e1f(ji,jj), e2f(ji,jj) )
61               ahm1(ji,jj) = za00 * zetmax
62               ahm2(ji,jj) = za00 * zefmax
63            END DO
64         END DO
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         ! Special case for ORCA R1, R2 and R4 configurations (overwrite the value of ahm1 ahm2)
74         ! ==============================================
75         IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) )   CALL ldf_dyn_c2d_orca( ld_print )
76         IF( cp_cfg == "orca" .AND.   jp_cfg == 1)                       CALL ldf_dyn_c2d_orca_R1( ld_print )
77
78         ! Control print
79         IF( lwp .AND. ld_print ) THEN
80            WRITE(numout,*)
81            WRITE(numout,*) 'inildf: 2D ahm1 array'
82            CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
83            WRITE(numout,*)
84            WRITE(numout,*) 'inildf: 2D ahm2 array'
85            CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
86         ENDIF
87      ENDIF
88
89      ! biharmonic operator (ahm3, ahm4) : at U- and V-points (used for bilaplacian operator
90      ! =================================                      whatever its orientation is)
91      IF( ln_dynldf_bilap ) THEN
92         ! (USER: modify ahm3 and ahm4 following your desiderata)
93         ! Here: ahm is proportional to the cube of the maximum of the gridspacing
94         !       in the to horizontal direction
95
96         zd_max = MAX( MAXVAL( e1u(:,:) ), MAXVAL( e2u(:,:) ) )
97         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
98
99         IF(lwp) WRITE(numout,*) '              bi-laplacian operator: ahm proportional to e1**3 '
100         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
101
102         za00 = ahm0_blp / ( zd_max * zd_max * zd_max )
103         DO jj = 1, jpj
104            DO ji = 1, jpi
105               zeumax = MAX( e1u(ji,jj), e2u(ji,jj) )
106               zevmax = MAX( e1v(ji,jj), e2v(ji,jj) )
107               ahm3(ji,jj) = za00 * zeumax * zeumax * zeumax
108               ahm4(ji,jj) = za00 * zevmax * zevmax * zevmax
109            END DO
110         END DO
111
112         ! Control print
113         IF( lwp .AND. ld_print ) THEN
114            WRITE(numout,*)
115            WRITE(numout,*) 'inildf: ahm3 array'
116            CALL prihre(ahm3,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
117            WRITE(numout,*)
118            WRITE(numout,*) 'inildf: ahm4 array'
119            CALL prihre(ahm4,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
120         ENDIF
121      ENDIF
122      !
123   END SUBROUTINE ldf_dyn_c2d
124
125
126   SUBROUTINE ldf_dyn_c2d_orca( ld_print )
127      !!----------------------------------------------------------------------
128      !!                 ***  ROUTINE ldf_dyn_c2d  ***
129      !!
130      !!                   **** W A R N I N G ****
131      !!
132      !!                ORCA R2 and R4 configurations
133      !!                 
134      !!                   **** W A R N I N G ****
135      !!                 
136      !! ** Purpose :   initializations of the lateral viscosity for orca R2
137      !!
138      !! ** Method  :   blah blah blah...
139      !!
140      !!----------------------------------------------------------------------
141      USE ldftra_oce, ONLY:   aht0
142      !
143      LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout
144      !
145      INTEGER  ::   ji, jj, jn   ! dummy loop indices
146      INTEGER  ::   inum, iim, ijm            ! local integers
147      INTEGER  ::   ifreq, il1, il2, ij, ii
148      INTEGER  ::   ijpt0,ijpt1
149      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk
150      CHARACTER (len=15) ::   clexp
151      INTEGER, POINTER, DIMENSION(:,:)  :: icof
152      INTEGER, POINTER, DIMENSION(:,:)  :: idata
153      !!----------------------------------------------------------------------
154      !                               
155      CALL wrk_alloc( jpi   , jpj   , icof  )
156      CALL wrk_alloc( jpidta, jpjdta, idata )
157      !
158      IF(lwp) WRITE(numout,*)
159      IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
160      IF(lwp) WRITE(numout,*) '~~~~~~  --'
161      IF(lwp) WRITE(numout,*) '        orca ocean configuration'
162
163      IF( cp_cfg == "orca" .AND. cp_cfz == "antarctic" ) THEN
164!
165! 1.2 Modify ahm
166! --------------
167         IF(lwp)WRITE(numout,*) ' inildf: Antarctic ocean'
168         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
169         IF(lwp)WRITE(numout,*) '         north boundary increase'
170
171         ahm1(:,:) = ahm0
172         ahm2(:,:) = ahm0
173
174         ijpt0=max(1,min(49 -njmpp+1,jpj))
175         ijpt1=max(0,min(49-njmpp+1,jpj-1))
176         DO jj=ijpt0,ijpt1
177            ahm2(:,jj)=ahm0*2.
178            ahm1(:,jj)=ahm0*2.
179         END DO
180         ijpt0=max(1,min(48 -njmpp+1,jpj))
181         ijpt1=max(0,min(48-njmpp+1,jpj-1))
182         DO jj=ijpt0,ijpt1
183            ahm2(:,jj)=ahm0*1.9
184            ahm1(:,jj)=ahm0*1.75
185         END DO
186         ijpt0=max(1,min(47 -njmpp+1,jpj))
187         ijpt1=max(0,min(47-njmpp+1,jpj-1))
188         DO jj=ijpt0,ijpt1
189            ahm2(:,jj)=ahm0*1.5
190            ahm1(:,jj)=ahm0*1.25
191         END DO
192         ijpt0=max(1,min(46 -njmpp+1,jpj))
193         ijpt1=max(0,min(46-njmpp+1,jpj-1))
194         DO jj=ijpt0,ijpt1
195            ahm2(:,jj)=ahm0*1.1
196         END DO
197
198      ELSE IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
199! 1.2 Modify ahm
200! --------------
201         IF(lwp)WRITE(numout,*) ' inildf: Arctic ocean'
202         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
203         IF(lwp)WRITE(numout,*) '         south and west boundary increase'
204
205
206         ahm1(:,:) = ahm0
207         ahm2(:,:) = ahm0
208
209         ijpt0=max(1,min(98-jpjzoom+1-njmpp+1,jpj))
210         ijpt1=max(0,min(98-jpjzoom+1-njmpp+1,jpj-1))
211         DO jj=ijpt0,ijpt1
212            ahm2(:,jj)=ahm0*2.
213            ahm1(:,jj)=ahm0*2.
214         END DO
215         ijpt0=max(1,min(99-jpjzoom+1-njmpp+1,jpj))
216         ijpt1=max(0,min(99-jpjzoom+1-njmpp+1,jpj-1))
217         DO jj=ijpt0,ijpt1
218            ahm2(:,jj)=ahm0*1.9
219            ahm1(:,jj)=ahm0*1.75
220         END DO
221         ijpt0=max(1,min(100-jpjzoom+1-njmpp+1,jpj))
222         ijpt1=max(0,min(100-jpjzoom+1-njmpp+1,jpj-1))
223         DO jj=ijpt0,ijpt1
224            ahm2(:,jj)=ahm0*1.5
225            ahm1(:,jj)=ahm0*1.25
226         END DO
227         ijpt0=max(1,min(101-jpjzoom+1-njmpp+1,jpj))
228         ijpt1=max(0,min(101-jpjzoom+1-njmpp+1,jpj-1))
229         DO jj=ijpt0,ijpt1
230            ahm2(:,jj)=ahm0*1.1
231         END DO
232      ELSE
233         ! Read 2d integer array to specify western boundary increase in the
234         ! ===================== equatorial strip (20N-20S) defined at t-points
235         
236         CALL ctl_opn( inum, 'ahmcoef', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
237         READ(inum,9101) clexp, iim, ijm
238         READ(inum,'(/)')
239         ifreq = 40
240         il1 = 1
241         DO jn = 1, jpidta/ifreq+1
242            READ(inum,'(/)')
243            il2 = MIN( jpidta, il1+ifreq-1 )
244            READ(inum,9201) ( ii, ji = il1, il2, 5 )
245            READ(inum,'(/)')
246            DO jj = jpjdta, 1, -1
247               READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
248            END DO
249            il1 = il1 + ifreq
250         END DO
251
252         DO jj = 1, nlcj
253            DO ji = 1, nlci
254               icof(ji,jj) = idata( mig(ji), mjg(jj) )
255            END DO
256         END DO
257         DO jj = nlcj+1, jpj
258            DO ji = 1, nlci
259               icof(ji,jj) = icof(ji,nlcj)
260            END DO
261         END DO
262         DO jj = 1, jpj
263            DO ji = nlci+1, jpi
264               icof(ji,jj) = icof(nlci,jj)
265            END DO
266         END DO
267
2689101     FORMAT(1x,a15,2i8)
2699201     FORMAT(3x,13(i3,12x))
2709202     FORMAT(i3,41i3)
271
272
273         ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator)
274         ! =================
275         ! define ahm1 and ahm2 at the right grid point position
276         ! (USER: modify ahm1 and ahm2 following your desiderata)
277
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         zahmeq = aht0
286
287         DO jj = 1, jpj
288            DO ji = 1, jpi
289               IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
290                  ahm2(ji,jj) =  ahm0
291               ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
292                  ahm2(ji,jj) =  zahmeq
293               ELSE
294                  ahm2(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   &
295                     * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
296               ENDIF
297               IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
298                  ahm1(ji,jj) =  ahm0
299               ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
300                  ahm1(ji,jj) =  zahmeq
301               ELSE
302                  ahm1(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   &
303                     * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
304               ENDIF
305            END DO
306         END DO
307
308         ! increase along western boundaries of equatorial strip
309         ! t-point
310         DO jj = 1, jpjm1
311            DO ji = 1, jpim1
312               zcoft = FLOAT( icof(ji,jj) ) / 100.
313               ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
314            END DO
315         END DO
316         ! f-point
317         icof(:,:) = icof(:,:) * tmask(:,:,1)
318         DO jj = 1, jpjm1
319            DO ji = 1, jpim1   ! NO vector opt.
320               zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
321               IF( zmsk == 0. ) THEN
322                  zcoff = 1.
323               ELSE
324                  zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   &
325                     / (zmsk * 100.)
326               ENDIF
327               ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
328            END DO
329         END DO
330      ENDIF
331     
332      ! Lateral boundary conditions on ( ahm1, ahm2 )
333      !                                ==============
334      CALL lbc_lnk( ahm1, 'T', 1. )   ! T-point, unchanged sign
335      CALL lbc_lnk( ahm2, 'F', 1. )   ! F-point, unchanged sign
336
337      ! Control print
338      IF( lwp .AND. ld_print ) THEN
339         WRITE(numout,*)
340         WRITE(numout,*) 'inildf: 2D ahm1 array'
341         CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
342         WRITE(numout,*)
343         WRITE(numout,*) 'inildf: 2D ahm2 array'
344         CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
345      ENDIF
346      !
347      CALL wrk_dealloc( jpi   , jpj   , icof  )
348      CALL wrk_dealloc( jpidta, jpjdta, idata )
349      !
350   END SUBROUTINE ldf_dyn_c2d_orca
351
352
353   SUBROUTINE ldf_dyn_c2d_orca_R1( ld_print )
354      !!----------------------------------------------------------------------
355      !!                 ***  ROUTINE ldf_dyn_c2d  ***
356      !!
357      !!                   **** W A R N I N G ****
358      !!
359      !!                ORCA R1 configuration
360      !!                 
361      !!                   **** W A R N I N G ****
362      !!                 
363      !! ** Purpose :   initializations of the lateral viscosity for orca R1
364      !!
365      !! ** Method  :   blah blah blah...
366      !!
367      !!----------------------------------------------------------------------
368      USE ldftra_oce, ONLY:   aht0
369      !
370      LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout
371      !
372      INTEGER ::   ji, jj, jn      ! dummy loop indices
373      INTEGER ::   inum            ! temporary logical unit
374      INTEGER ::   iim, ijm
375      INTEGER ::   ifreq, il1, il2, ij, ii
376      INTEGER ::   ijpt0,ijpt1
377      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk, zam20s
378      CHARACTER (len=15) ::   clexp
379      INTEGER, POINTER, DIMENSION(:,:)  :: icof
380      INTEGER, POINTER, DIMENSION(:,:)  :: idata
381      !!----------------------------------------------------------------------
382      !                               
383      CALL wrk_alloc( jpi   , jpj   , icof  )
384      CALL wrk_alloc( jpidta, jpjdta, idata )
385      !                               
386
387      IF(lwp) WRITE(numout,*)
388      IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
389      IF(lwp) WRITE(numout,*) '~~~~~~  --'
390      IF(lwp) WRITE(numout,*) '        orca_r1 configuration'
391
392      IF( cp_cfg == "orca" .AND. cp_cfz == "antarctic" ) THEN
393!
394! 1.2 Modify ahm
395! --------------
396         IF(lwp)WRITE(numout,*) ' inildf: Antarctic ocean'
397         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
398         IF(lwp)WRITE(numout,*) '         north boundary increase'
399
400         ahm1(:,:) = ahm0
401         ahm2(:,:) = ahm0
402
403         ijpt0=max(1,min(49 -njmpp+1,jpj))
404         ijpt1=max(0,min(49-njmpp+1,jpj-1))
405         DO jj=ijpt0,ijpt1
406            ahm2(:,jj)=ahm0*2.
407            ahm1(:,jj)=ahm0*2.
408         END DO
409         ijpt0=max(1,min(48 -njmpp+1,jpj))
410         ijpt1=max(0,min(48-njmpp+1,jpj-1))
411         DO jj=ijpt0,ijpt1
412            ahm2(:,jj)=ahm0*1.9
413            ahm1(:,jj)=ahm0*1.75
414         END DO
415         ijpt0=max(1,min(47 -njmpp+1,jpj))
416         ijpt1=max(0,min(47-njmpp+1,jpj-1))
417         DO jj=ijpt0,ijpt1
418            ahm2(:,jj)=ahm0*1.5
419            ahm1(:,jj)=ahm0*1.25
420         END DO
421         ijpt0=max(1,min(46 -njmpp+1,jpj))
422         ijpt1=max(0,min(46-njmpp+1,jpj-1))
423         DO jj=ijpt0,ijpt1
424            ahm2(:,jj)=ahm0*1.1
425         END DO
426
427      ELSE IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
428! 1.2 Modify ahm
429! --------------
430         IF(lwp)WRITE(numout,*) ' inildf: Arctic ocean'
431         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
432         IF(lwp)WRITE(numout,*) '         south and west boundary increase'
433
434
435         ahm1(:,:) = ahm0
436         ahm2(:,:) = ahm0
437
438         ijpt0=max(1,min(98-jpjzoom+1-njmpp+1,jpj))
439         ijpt1=max(0,min(98-jpjzoom+1-njmpp+1,jpj-1))
440         DO jj=ijpt0,ijpt1
441            ahm2(:,jj)=ahm0*2.
442            ahm1(:,jj)=ahm0*2.
443         END DO
444         ijpt0=max(1,min(99-jpjzoom+1-njmpp+1,jpj))
445         ijpt1=max(0,min(99-jpjzoom+1-njmpp+1,jpj-1))
446         DO jj=ijpt0,ijpt1
447            ahm2(:,jj)=ahm0*1.9
448            ahm1(:,jj)=ahm0*1.75
449         END DO
450         ijpt0=max(1,min(100-jpjzoom+1-njmpp+1,jpj))
451         ijpt1=max(0,min(100-jpjzoom+1-njmpp+1,jpj-1))
452         DO jj=ijpt0,ijpt1
453            ahm2(:,jj)=ahm0*1.5
454            ahm1(:,jj)=ahm0*1.25
455         END DO
456         ijpt0=max(1,min(101-jpjzoom+1-njmpp+1,jpj))
457         ijpt1=max(0,min(101-jpjzoom+1-njmpp+1,jpj-1))
458         DO jj=ijpt0,ijpt1
459            ahm2(:,jj)=ahm0*1.1
460         END DO
461      ELSE
462         
463         ! Read 2d integer array to specify western boundary increase in the
464         ! ===================== equatorial strip (20N-20S) defined at t-points
465         
466         CALL ctl_opn( inum, 'ahmcoef', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
467            &           1, numout, lwp )
468         REWIND inum
469         READ(inum,9101) clexp, iim, ijm
470         READ(inum,'(/)')
471         ifreq = 40
472         il1 = 1
473         DO jn = 1, jpidta/ifreq+1
474            READ(inum,'(/)')
475            il2 = MIN( jpidta, il1+ifreq-1 )
476            READ(inum,9201) ( ii, ji = il1, il2, 5 )
477            READ(inum,'(/)')
478            DO jj = jpjdta, 1, -1
479               READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
480            END DO
481            il1 = il1 + ifreq
482         END DO
483
484         DO jj = 1, nlcj
485            DO ji = 1, nlci
486               icof(ji,jj) = idata( mig(ji), mjg(jj) )
487            END DO
488         END DO
489         DO jj = nlcj+1, jpj
490            DO ji = 1, nlci
491               icof(ji,jj) = icof(ji,nlcj)
492            END DO
493         END DO
494         DO jj = 1, jpj
495            DO ji = nlci+1, jpi
496               icof(ji,jj) = icof(nlci,jj)
497            END DO
498         END DO
499
5009101     FORMAT(1x,a15,2i8)
5019201     FORMAT(3x,13(i3,12x))
5029202     FORMAT(i3,41i3)
503
504
505         ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator)
506         ! =================
507         ! define ahm1 and ahm2 at the right grid point position
508         ! (USER: modify ahm1 and ahm2 following your desiderata)
509
510
511         ! Decrease ahm to zahmeq m2/s in the tropics
512         ! (from 90   to 20   degrees: ahm = scaled by local metrics
513         !  from 20   to  2.5 degrees: ahm = decrease in (1-cos)/2
514         !  from  2.5 to  0   degrees: ahm = constant
515         ! symmetric in the south hemisphere)
516
517         zahmeq = aht0
518         zam20s = ahm0*COS( rad * 20. )
519
520         DO jj = 1, jpj
521            DO ji = 1, jpi
522               IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
523                  !              leave as set in ldf_dyn_c2d
524               ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
525                  ahm2(ji,jj) =  zahmeq
526               ELSE
527                  ahm2(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   &
528                     * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
529               ENDIF
530               IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
531                  !             leave as set in ldf_dyn_c2d
532               ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
533                  ahm1(ji,jj) =  zahmeq
534               ELSE
535                  ahm1(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   &
536                     * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
537               ENDIF
538            END DO
539         END DO
540
541         ! increase along western boundaries of equatorial strip
542         ! t-point
543         DO jj = 1, jpjm1
544            DO ji = 1, jpim1
545               IF( ABS( gphit(ji,jj) ) < 20. ) THEN
546                  zcoft = FLOAT( icof(ji,jj) ) / 100.
547                  ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
548               ENDIF
549            END DO
550         END DO
551         ! f-point
552         icof(:,:) = icof(:,:) * tmask(:,:,1)
553         DO jj = 1, jpjm1
554            DO ji = 1, jpim1
555               IF( ABS( gphif(ji,jj) ) < 20. ) THEN
556                  zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
557                  IF( zmsk == 0. ) THEN
558                     zcoff = 1.
559                  ELSE
560                     zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   &
561                        / (zmsk * 100.)
562                  ENDIF
563                  ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
564               ENDIF
565            END DO
566         END DO
567      ENDIF
568     
569      ! Lateral boundary conditions on ( ahm1, ahm2 )
570      !                                ==============
571      CALL lbc_lnk( ahm1, 'T', 1. )   ! T-point, unchanged sign
572      CALL lbc_lnk( ahm2, 'F', 1. )   ! F-point, unchanged sign
573
574      ! Control print
575      IF( lwp .AND. ld_print ) THEN
576         WRITE(numout,*)
577         WRITE(numout,*) 'inildf: 2D ahm1 array'
578         CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
579         WRITE(numout,*)
580         WRITE(numout,*) 'inildf: 2D ahm2 array'
581         CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
582      ENDIF
583      !
584      CALL wrk_dealloc( jpi   , jpj   , icof  )
585      CALL wrk_dealloc( jpidta, jpjdta, idata )
586      !
587   END SUBROUTINE ldf_dyn_c2d_orca_R1
Note: See TracBrowser for help on using the repository browser.