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

source: branches/2013/dev_r3853_CNRS9_ConfSetting/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90 @ 3989

Last change on this file since 3989 was 3989, checked in by clevy, 11 years ago

Configuration setting/Step3 and doc, see ticket:#1074
  • Property svn:keywords set to Id
File size: 21.8 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      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk
149      CHARACTER (len=15) ::   clexp
150      INTEGER, POINTER, DIMENSION(:,:)  :: icof
151      INTEGER, POINTER, DIMENSION(:,:)  :: idata
152      !!----------------------------------------------------------------------
153      !                               
154      CALL wrk_alloc( jpi   , jpj   , icof  )
155      CALL wrk_alloc( jpidta, jpjdta, idata )
156      !
157      IF(lwp) WRITE(numout,*)
158      IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
159      IF(lwp) WRITE(numout,*) '~~~~~~  --'
160      IF(lwp) WRITE(numout,*) '        orca ocean configuration'
161
162      IF( cp_cfg == "orca" .AND. cp_cfz == "antarctic" ) THEN
163!
164! 1.2 Modify ahm
165! --------------
166         IF(lwp)WRITE(numout,*) ' inildf: Antarctic ocean'
167         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
168         IF(lwp)WRITE(numout,*) '         north boundary increase'
169
170         ahm1(:,:) = ahm0
171         ahm2(:,:) = ahm0
172
173         ijpt0=max(1,min(49 -njmpp+1,jpj))
174         ijpt1=max(0,min(49-njmpp+1,jpj-1))
175         DO jj=ijpt0,ijpt1
176            ahm2(:,jj)=ahm0*2.
177            ahm1(:,jj)=ahm0*2.
178         END DO
179         ijpt0=max(1,min(48 -njmpp+1,jpj))
180         ijpt1=max(0,min(48-njmpp+1,jpj-1))
181         DO jj=ijpt0,ijpt1
182            ahm2(:,jj)=ahm0*1.9
183            ahm1(:,jj)=ahm0*1.75
184         END DO
185         ijpt0=max(1,min(47 -njmpp+1,jpj))
186         ijpt1=max(0,min(47-njmpp+1,jpj-1))
187         DO jj=ijpt0,ijpt1
188            ahm2(:,jj)=ahm0*1.5
189            ahm1(:,jj)=ahm0*1.25
190         END DO
191         ijpt0=max(1,min(46 -njmpp+1,jpj))
192         ijpt1=max(0,min(46-njmpp+1,jpj-1))
193         DO jj=ijpt0,ijpt1
194            ahm2(:,jj)=ahm0*1.1
195         END DO
196
197      ELSE IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
198! 1.2 Modify ahm
199! --------------
200         IF(lwp)WRITE(numout,*) ' inildf: Arctic ocean'
201         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
202         IF(lwp)WRITE(numout,*) '         south and west boundary increase'
203
204
205         ahm1(:,:) = ahm0
206         ahm2(:,:) = ahm0
207
208         ijpt0=max(1,min(98-jpjzoom+1-njmpp+1,jpj))
209         ijpt1=max(0,min(98-jpjzoom+1-njmpp+1,jpj-1))
210         DO jj=ijpt0,ijpt1
211            ahm2(:,jj)=ahm0*2.
212            ahm1(:,jj)=ahm0*2.
213         END DO
214         ijpt0=max(1,min(99-jpjzoom+1-njmpp+1,jpj))
215         ijpt1=max(0,min(99-jpjzoom+1-njmpp+1,jpj-1))
216         DO jj=ijpt0,ijpt1
217            ahm2(:,jj)=ahm0*1.9
218            ahm1(:,jj)=ahm0*1.75
219         END DO
220         ijpt0=max(1,min(100-jpjzoom+1-njmpp+1,jpj))
221         ijpt1=max(0,min(100-jpjzoom+1-njmpp+1,jpj-1))
222         DO jj=ijpt0,ijpt1
223            ahm2(:,jj)=ahm0*1.5
224            ahm1(:,jj)=ahm0*1.25
225         END DO
226         ijpt0=max(1,min(101-jpjzoom+1-njmpp+1,jpj))
227         ijpt1=max(0,min(101-jpjzoom+1-njmpp+1,jpj-1))
228         DO jj=ijpt0,ijpt1
229            ahm2(:,jj)=ahm0*1.1
230         END DO
231      ELSE
232         ! Read 2d integer array to specify western boundary increase in the
233         ! ===================== equatorial strip (20N-20S) defined at t-points
234         
235         CALL ctl_opn( inum, 'ahmcoef', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
236         READ(inum,9101) clexp, iim, ijm
237         READ(inum,'(/)')
238         ifreq = 40
239         il1 = 1
240         DO jn = 1, jpidta/ifreq+1
241            READ(inum,'(/)')
242            il2 = MIN( jpidta, il1+ifreq-1 )
243            READ(inum,9201) ( ii, ji = il1, il2, 5 )
244            READ(inum,'(/)')
245            DO jj = jpjdta, 1, -1
246               READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
247            END DO
248            il1 = il1 + ifreq
249         END DO
250
251         DO jj = 1, nlcj
252            DO ji = 1, nlci
253               icof(ji,jj) = idata( mig(ji), mjg(jj) )
254            END DO
255         END DO
256         DO jj = nlcj+1, jpj
257            DO ji = 1, nlci
258               icof(ji,jj) = icof(ji,nlcj)
259            END DO
260         END DO
261         DO jj = 1, jpj
262            DO ji = nlci+1, jpi
263               icof(ji,jj) = icof(nlci,jj)
264            END DO
265         END DO
266
2679101     FORMAT(1x,a15,2i8)
2689201     FORMAT(3x,13(i3,12x))
2699202     FORMAT(i3,41i3)
270
271
272         ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator)
273         ! =================
274         ! define ahm1 and ahm2 at the right grid point position
275         ! (USER: modify ahm1 and ahm2 following your desiderata)
276
277
278         ! Decrease ahm to zahmeq m2/s in the tropics
279         ! (from 90 to 20 degre: ahm = constant
280         ! from 20 to  2.5 degre: ahm = decrease in (1-cos)/2
281         ! from  2.5 to  0 degre: ahm = constant
282         ! symmetric in the south hemisphere)
283
284         zahmeq = aht0
285
286         DO jj = 1, jpj
287            DO ji = 1, jpi
288               IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
289                  ahm2(ji,jj) =  ahm0
290               ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
291                  ahm2(ji,jj) =  zahmeq
292               ELSE
293                  ahm2(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   &
294                     * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
295               ENDIF
296               IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
297                  ahm1(ji,jj) =  ahm0
298               ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
299                  ahm1(ji,jj) =  zahmeq
300               ELSE
301                  ahm1(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   &
302                     * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
303               ENDIF
304            END DO
305         END DO
306
307         ! increase along western boundaries of equatorial strip
308         ! t-point
309         DO jj = 1, jpjm1
310            DO ji = 1, jpim1
311               zcoft = FLOAT( icof(ji,jj) ) / 100.
312               ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
313            END DO
314         END DO
315         ! f-point
316         icof(:,:) = icof(:,:) * tmask(:,:,1)
317         DO jj = 1, jpjm1
318            DO ji = 1, jpim1   ! NO vector opt.
319               zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
320               IF( zmsk == 0. ) THEN
321                  zcoff = 1.
322               ELSE
323                  zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   &
324                     / (zmsk * 100.)
325               ENDIF
326               ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
327            END DO
328         END DO
329      ENDIF
330     
331      ! Lateral boundary conditions on ( ahm1, ahm2 )
332      !                                ==============
333      CALL lbc_lnk( ahm1, 'T', 1. )   ! T-point, unchanged sign
334      CALL lbc_lnk( ahm2, 'F', 1. )   ! F-point, unchanged sign
335
336      ! Control print
337      IF( lwp .AND. ld_print ) THEN
338         WRITE(numout,*)
339         WRITE(numout,*) 'inildf: 2D ahm1 array'
340         CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
341         WRITE(numout,*)
342         WRITE(numout,*) 'inildf: 2D ahm2 array'
343         CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
344      ENDIF
345      !
346      CALL wrk_dealloc( jpi   , jpj   , icof  )
347      CALL wrk_dealloc( jpidta, jpjdta, idata )
348      !
349   END SUBROUTINE ldf_dyn_c2d_orca
350
351
352   SUBROUTINE ldf_dyn_c2d_orca_R1( ld_print )
353      !!----------------------------------------------------------------------
354      !!                 ***  ROUTINE ldf_dyn_c2d  ***
355      !!
356      !!                   **** W A R N I N G ****
357      !!
358      !!                ORCA R1 configuration
359      !!                 
360      !!                   **** W A R N I N G ****
361      !!                 
362      !! ** Purpose :   initializations of the lateral viscosity for orca R1
363      !!
364      !! ** Method  :   blah blah blah...
365      !!
366      !!----------------------------------------------------------------------
367      USE ldftra_oce, ONLY:   aht0
368      !
369      LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout
370      !
371      INTEGER ::   ji, jj, jn      ! dummy loop indices
372      INTEGER ::   inum            ! temporary logical unit
373      INTEGER ::   iim, ijm
374      INTEGER ::   ifreq, il1, il2, ij, ii
375      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk, zam20s
376      CHARACTER (len=15) ::   clexp
377      INTEGER, POINTER, DIMENSION(:,:)  :: icof
378      INTEGER, POINTER, DIMENSION(:,:)  :: idata
379      !!----------------------------------------------------------------------
380      !                               
381      CALL wrk_alloc( jpi   , jpj   , icof  )
382      CALL wrk_alloc( jpidta, jpjdta, idata )
383      !                               
384
385      IF(lwp) WRITE(numout,*)
386      IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
387      IF(lwp) WRITE(numout,*) '~~~~~~  --'
388      IF(lwp) WRITE(numout,*) '        orca_r1 configuration'
389
390      IF( cp_cfg == "orca" .AND. cp_cfz == "antarctic" ) THEN
391!
392! 1.2 Modify ahm
393! --------------
394         IF(lwp)WRITE(numout,*) ' inildf: Antarctic ocean'
395         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
396         IF(lwp)WRITE(numout,*) '         north boundary increase'
397
398         ahm1(:,:) = ahm0
399         ahm2(:,:) = ahm0
400
401         ijpt0=max(1,min(49 -njmpp+1,jpj))
402         ijpt1=max(0,min(49-njmpp+1,jpj-1))
403         DO jj=ijpt0,ijpt1
404            ahm2(:,jj)=ahm0*2.
405            ahm1(:,jj)=ahm0*2.
406         END DO
407         ijpt0=max(1,min(48 -njmpp+1,jpj))
408         ijpt1=max(0,min(48-njmpp+1,jpj-1))
409         DO jj=ijpt0,ijpt1
410            ahm2(:,jj)=ahm0*1.9
411            ahm1(:,jj)=ahm0*1.75
412         END DO
413         ijpt0=max(1,min(47 -njmpp+1,jpj))
414         ijpt1=max(0,min(47-njmpp+1,jpj-1))
415         DO jj=ijpt0,ijpt1
416            ahm2(:,jj)=ahm0*1.5
417            ahm1(:,jj)=ahm0*1.25
418         END DO
419         ijpt0=max(1,min(46 -njmpp+1,jpj))
420         ijpt1=max(0,min(46-njmpp+1,jpj-1))
421         DO jj=ijpt0,ijpt1
422            ahm2(:,jj)=ahm0*1.1
423         END DO
424
425      ELSE IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
426! 1.2 Modify ahm
427! --------------
428         IF(lwp)WRITE(numout,*) ' inildf: Arctic ocean'
429         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
430         IF(lwp)WRITE(numout,*) '         south and west boundary increase'
431
432
433         ahm1(:,:) = ahm0
434         ahm2(:,:) = ahm0
435
436         ijpt0=max(1,min(98-jpjzoom+1-njmpp+1,jpj))
437         ijpt1=max(0,min(98-jpjzoom+1-njmpp+1,jpj-1))
438         DO jj=ijpt0,ijpt1
439            ahm2(:,jj)=ahm0*2.
440            ahm1(:,jj)=ahm0*2.
441         END DO
442         ijpt0=max(1,min(99-jpjzoom+1-njmpp+1,jpj))
443         ijpt1=max(0,min(99-jpjzoom+1-njmpp+1,jpj-1))
444         DO jj=ijpt0,ijpt1
445            ahm2(:,jj)=ahm0*1.9
446            ahm1(:,jj)=ahm0*1.75
447         END DO
448         ijpt0=max(1,min(100-jpjzoom+1-njmpp+1,jpj))
449         ijpt1=max(0,min(100-jpjzoom+1-njmpp+1,jpj-1))
450         DO jj=ijpt0,ijpt1
451            ahm2(:,jj)=ahm0*1.5
452            ahm1(:,jj)=ahm0*1.25
453         END DO
454         ijpt0=max(1,min(101-jpjzoom+1-njmpp+1,jpj))
455         ijpt1=max(0,min(101-jpjzoom+1-njmpp+1,jpj-1))
456         DO jj=ijpt0,ijpt1
457            ahm2(:,jj)=ahm0*1.1
458         END DO
459      ELSE
460         
461         ! Read 2d integer array to specify western boundary increase in the
462         ! ===================== equatorial strip (20N-20S) defined at t-points
463         
464         CALL ctl_opn( inum, 'ahmcoef', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
465            &           1, numout, lwp )
466         REWIND inum
467         READ(inum,9101) clexp, iim, ijm
468         READ(inum,'(/)')
469         ifreq = 40
470         il1 = 1
471         DO jn = 1, jpidta/ifreq+1
472            READ(inum,'(/)')
473            il2 = MIN( jpidta, il1+ifreq-1 )
474            READ(inum,9201) ( ii, ji = il1, il2, 5 )
475            READ(inum,'(/)')
476            DO jj = jpjdta, 1, -1
477               READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
478            END DO
479            il1 = il1 + ifreq
480         END DO
481
482         DO jj = 1, nlcj
483            DO ji = 1, nlci
484               icof(ji,jj) = idata( mig(ji), mjg(jj) )
485            END DO
486         END DO
487         DO jj = nlcj+1, jpj
488            DO ji = 1, nlci
489               icof(ji,jj) = icof(ji,nlcj)
490            END DO
491         END DO
492         DO jj = 1, jpj
493            DO ji = nlci+1, jpi
494               icof(ji,jj) = icof(nlci,jj)
495            END DO
496         END DO
497
4989101     FORMAT(1x,a15,2i8)
4999201     FORMAT(3x,13(i3,12x))
5009202     FORMAT(i3,41i3)
501
502
503         ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator)
504         ! =================
505         ! define ahm1 and ahm2 at the right grid point position
506         ! (USER: modify ahm1 and ahm2 following your desiderata)
507
508
509         ! Decrease ahm to zahmeq m2/s in the tropics
510         ! (from 90   to 20   degrees: ahm = scaled by local metrics
511         !  from 20   to  2.5 degrees: ahm = decrease in (1-cos)/2
512         !  from  2.5 to  0   degrees: ahm = constant
513         ! symmetric in the south hemisphere)
514
515         zahmeq = aht0
516         zam20s = ahm0*COS( rad * 20. )
517
518         DO jj = 1, jpj
519            DO ji = 1, jpi
520               IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
521                  !              leave as set in ldf_dyn_c2d
522               ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
523                  ahm2(ji,jj) =  zahmeq
524               ELSE
525                  ahm2(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   &
526                     * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
527               ENDIF
528               IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
529                  !             leave as set in ldf_dyn_c2d
530               ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
531                  ahm1(ji,jj) =  zahmeq
532               ELSE
533                  ahm1(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   &
534                     * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
535               ENDIF
536            END DO
537         END DO
538
539         ! increase along western boundaries of equatorial strip
540         ! t-point
541         DO jj = 1, jpjm1
542            DO ji = 1, jpim1
543               IF( ABS( gphit(ji,jj) ) < 20. ) THEN
544                  zcoft = FLOAT( icof(ji,jj) ) / 100.
545                  ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
546               ENDIF
547            END DO
548         END DO
549         ! f-point
550         icof(:,:) = icof(:,:) * tmask(:,:,1)
551         DO jj = 1, jpjm1
552            DO ji = 1, jpim1
553               IF( ABS( gphif(ji,jj) ) < 20. ) THEN
554                  zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
555                  IF( zmsk == 0. ) THEN
556                     zcoff = 1.
557                  ELSE
558                     zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   &
559                        / (zmsk * 100.)
560                  ENDIF
561                  ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
562               ENDIF
563            END DO
564         END DO
565      ENDIF
566     
567      ! Lateral boundary conditions on ( ahm1, ahm2 )
568      !                                ==============
569      CALL lbc_lnk( ahm1, 'T', 1. )   ! T-point, unchanged sign
570      CALL lbc_lnk( ahm2, 'F', 1. )   ! F-point, unchanged sign
571
572      ! Control print
573      IF( lwp .AND. ld_print ) THEN
574         WRITE(numout,*)
575         WRITE(numout,*) 'inildf: 2D ahm1 array'
576         CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
577         WRITE(numout,*)
578         WRITE(numout,*) 'inildf: 2D ahm2 array'
579         CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
580      ENDIF
581      !
582      CALL wrk_dealloc( jpi   , jpj   , icof  )
583      CALL wrk_dealloc( jpidta, jpjdta, idata )
584      !
585   END SUBROUTINE ldf_dyn_c2d_orca_R1
Note: See TracBrowser for help on using the repository browser.