source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90 @ 6486

Last change on this file since 6486 was 6486, checked in by davestorkey, 5 years ago

Remove SVN keywords from UKMO/dev_r5518_GO6_package branch.

File size: 20.3 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      USE iom
143      !
144      LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout
145      !
146      INTEGER  ::   ji, jj, jn   ! dummy loop indices
147      INTEGER  ::   inum, iim, ijm            ! local integers
148      INTEGER  ::   ifreq, il1, il2, ij, ii
149      INTEGER  ::   ijpt0,ijpt1, ierror
150      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk
151      CHARACTER (len=15) ::   clexp
152      INTEGER,     POINTER, DIMENSION(:,:)  :: icof
153      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztemp2d  ! temporary array to read ahmcoef file
154      !!----------------------------------------------------------------------
155      !                               
156      CALL wrk_alloc( jpi   , jpj   , icof  )
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         ALLOCATE( ztemp2d(jpi,jpj) )
237         ztemp2d(:,:) = 0.
238         CALL iom_open ( 'ahmcoef.nc', inum )
239         CALL iom_get  ( inum, jpdom_data, 'icof', ztemp2d)
240         icof(:,:)  = NINT(ztemp2d(:,:))
241         CALL iom_close( inum )
242         DEALLOCATE(ztemp2d)
243
244         ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator)
245         ! =================
246         ! define ahm1 and ahm2 at the right grid point position
247         ! (USER: modify ahm1 and ahm2 following your desiderata)
248
249
250         ! Decrease ahm to zahmeq m2/s in the tropics
251         ! (from 90 to 20 degre: ahm = constant
252         ! from 20 to  2.5 degre: ahm = decrease in (1-cos)/2
253         ! from  2.5 to  0 degre: ahm = constant
254         ! symmetric in the south hemisphere)
255
256         zahmeq = aht0
257
258         DO jj = 1, jpj
259            DO ji = 1, jpi
260               IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
261                  ahm2(ji,jj) =  ahm0
262               ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
263                  ahm2(ji,jj) =  zahmeq
264               ELSE
265                  ahm2(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   &
266                     * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
267               ENDIF
268               IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
269                  ahm1(ji,jj) =  ahm0
270               ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
271                  ahm1(ji,jj) =  zahmeq
272               ELSE
273                  ahm1(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   &
274                     * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
275               ENDIF
276            END DO
277         END DO
278
279         ! increase along western boundaries of equatorial strip
280         ! t-point
281         DO jj = 1, jpjm1
282            DO ji = 1, jpim1
283               zcoft = FLOAT( icof(ji,jj) ) / 100.
284               ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
285            END DO
286         END DO
287         ! f-point
288         icof(:,:) = icof(:,:) * tmask(:,:,1)
289         DO jj = 1, jpjm1
290            DO ji = 1, jpim1   ! NO vector opt.
291               zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
292               IF( zmsk == 0. ) THEN
293                  zcoff = 1.
294               ELSE
295                  zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   &
296                     / (zmsk * 100.)
297               ENDIF
298               ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
299            END DO
300         END DO
301      ENDIF
302     
303      ! Lateral boundary conditions on ( ahm1, ahm2 )
304      !                                ==============
305      CALL lbc_lnk( ahm1, 'T', 1. )   ! T-point, unchanged sign
306      CALL lbc_lnk( ahm2, 'F', 1. )   ! F-point, unchanged sign
307
308      ! Control print
309      IF( lwp .AND. ld_print ) THEN
310         WRITE(numout,*)
311         WRITE(numout,*) 'inildf: 2D ahm1 array'
312         CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
313         WRITE(numout,*)
314         WRITE(numout,*) 'inildf: 2D ahm2 array'
315         CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
316      ENDIF
317      !
318      CALL wrk_dealloc( jpi   , jpj   , icof  )
319      !
320   END SUBROUTINE ldf_dyn_c2d_orca
321
322
323   SUBROUTINE ldf_dyn_c2d_orca_R1( ld_print )
324      !!----------------------------------------------------------------------
325      !!                 ***  ROUTINE ldf_dyn_c2d  ***
326      !!
327      !!                   **** W A R N I N G ****
328      !!
329      !!                ORCA R1 configuration
330      !!                 
331      !!                   **** W A R N I N G ****
332      !!                 
333      !! ** Purpose :   initializations of the lateral viscosity for orca R1
334      !!
335      !! ** Method  :   blah blah blah...
336      !!
337      !!----------------------------------------------------------------------
338      USE ldftra_oce, ONLY:   aht0
339      USE iom
340      !
341      LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout
342      !
343      INTEGER ::   ji, jj, jn      ! dummy loop indices
344      INTEGER ::   inum            ! temporary logical unit
345      INTEGER ::   iim, ijm
346      INTEGER ::   ifreq, il1, il2, ij, ii
347      INTEGER ::   ijpt0,ijpt1, ierror
348      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk, zam20s
349      CHARACTER (len=15) ::   clexp
350      INTEGER,     POINTER, DIMENSION(:,:)  :: icof
351      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztemp2d  ! temporary array to read ahmcoef file
352      !!----------------------------------------------------------------------
353      !                               
354      CALL wrk_alloc( jpi   , jpj   , icof  )
355      !                               
356      IF(lwp) WRITE(numout,*)
357      IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient'
358      IF(lwp) WRITE(numout,*) '~~~~~~  --'
359      IF(lwp) WRITE(numout,*) '        orca_r1 configuration'
360
361      IF( cp_cfg == "orca" .AND. cp_cfz == "antarctic" ) THEN
362!
363! 1.2 Modify ahm
364! --------------
365         IF(lwp)WRITE(numout,*) ' inildf: Antarctic ocean'
366         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
367         IF(lwp)WRITE(numout,*) '         north boundary increase'
368
369         ahm1(:,:) = ahm0
370         ahm2(:,:) = ahm0
371
372         ijpt0=max(1,min(49 -njmpp+1,jpj))
373         ijpt1=max(0,min(49-njmpp+1,jpj-1))
374         DO jj=ijpt0,ijpt1
375            ahm2(:,jj)=ahm0*2.
376            ahm1(:,jj)=ahm0*2.
377         END DO
378         ijpt0=max(1,min(48 -njmpp+1,jpj))
379         ijpt1=max(0,min(48-njmpp+1,jpj-1))
380         DO jj=ijpt0,ijpt1
381            ahm2(:,jj)=ahm0*1.9
382            ahm1(:,jj)=ahm0*1.75
383         END DO
384         ijpt0=max(1,min(47 -njmpp+1,jpj))
385         ijpt1=max(0,min(47-njmpp+1,jpj-1))
386         DO jj=ijpt0,ijpt1
387            ahm2(:,jj)=ahm0*1.5
388            ahm1(:,jj)=ahm0*1.25
389         END DO
390         ijpt0=max(1,min(46 -njmpp+1,jpj))
391         ijpt1=max(0,min(46-njmpp+1,jpj-1))
392         DO jj=ijpt0,ijpt1
393            ahm2(:,jj)=ahm0*1.1
394         END DO
395
396      ELSE IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
397! 1.2 Modify ahm
398! --------------
399         IF(lwp)WRITE(numout,*) ' inildf: Arctic ocean'
400         IF(lwp)WRITE(numout,*) '         no tropics, no reduction of ahm'
401         IF(lwp)WRITE(numout,*) '         south and west boundary increase'
402
403
404         ahm1(:,:) = ahm0
405         ahm2(:,:) = ahm0
406
407         ijpt0=max(1,min(98-jpjzoom+1-njmpp+1,jpj))
408         ijpt1=max(0,min(98-jpjzoom+1-njmpp+1,jpj-1))
409         DO jj=ijpt0,ijpt1
410            ahm2(:,jj)=ahm0*2.
411            ahm1(:,jj)=ahm0*2.
412         END DO
413         ijpt0=max(1,min(99-jpjzoom+1-njmpp+1,jpj))
414         ijpt1=max(0,min(99-jpjzoom+1-njmpp+1,jpj-1))
415         DO jj=ijpt0,ijpt1
416            ahm2(:,jj)=ahm0*1.9
417            ahm1(:,jj)=ahm0*1.75
418         END DO
419         ijpt0=max(1,min(100-jpjzoom+1-njmpp+1,jpj))
420         ijpt1=max(0,min(100-jpjzoom+1-njmpp+1,jpj-1))
421         DO jj=ijpt0,ijpt1
422            ahm2(:,jj)=ahm0*1.5
423            ahm1(:,jj)=ahm0*1.25
424         END DO
425         ijpt0=max(1,min(101-jpjzoom+1-njmpp+1,jpj))
426         ijpt1=max(0,min(101-jpjzoom+1-njmpp+1,jpj-1))
427         DO jj=ijpt0,ijpt1
428            ahm2(:,jj)=ahm0*1.1
429         END DO
430      ELSE
431         
432         ! Read 2d integer array to specify western boundary increase in the
433         ! ===================== equatorial strip (20N-20S) defined at t-points
434         ALLOCATE( ztemp2d(jpi,jpj) )
435         ztemp2d(:,:) = 0.
436         CALL iom_open ( 'ahmcoef.nc', inum )
437         CALL iom_get  ( inum, jpdom_data, 'icof', ztemp2d)
438         icof(:,:)  = NINT(ztemp2d(:,:))
439         CALL iom_close( inum )
440         DEALLOCATE(ztemp2d)
441
442         ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator)
443         ! =================
444         ! define ahm1 and ahm2 at the right grid point position
445         ! (USER: modify ahm1 and ahm2 following your desiderata)
446
447
448         ! Decrease ahm to zahmeq m2/s in the tropics
449         ! (from 90   to 20   degrees: ahm = scaled by local metrics
450         !  from 20   to  2.5 degrees: ahm = decrease in (1-cos)/2
451         !  from  2.5 to  0   degrees: ahm = constant
452         ! symmetric in the south hemisphere)
453
454         zahmeq = aht0
455         zam20s = ahm0*COS( rad * 20. )
456
457         DO jj = 1, jpj
458            DO ji = 1, jpi
459               IF( ABS( gphif(ji,jj) ) >= 20. ) THEN
460                  !              leave as set in ldf_dyn_c2d
461               ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN
462                  ahm2(ji,jj) =  zahmeq
463               ELSE
464                  ahm2(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   &
465                     * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) )
466               ENDIF
467               IF( ABS( gphit(ji,jj) ) >= 20. ) THEN
468                  !             leave as set in ldf_dyn_c2d
469               ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN
470                  ahm1(ji,jj) =  zahmeq
471               ELSE
472                  ahm1(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   &
473                     * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) )
474               ENDIF
475            END DO
476         END DO
477
478         ! increase along western boundaries of equatorial strip
479         ! t-point
480         DO jj = 1, jpjm1
481            DO ji = 1, jpim1
482               IF( ABS( gphit(ji,jj) ) < 20. ) THEN
483                  zcoft = FLOAT( icof(ji,jj) ) / 100.
484                  ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)
485               ENDIF
486            END DO
487         END DO
488         ! f-point
489         icof(:,:) = icof(:,:) * tmask(:,:,1)
490         DO jj = 1, jpjm1
491            DO ji = 1, jpim1
492               IF( ABS( gphif(ji,jj) ) < 20. ) THEN
493                  zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
494                  IF( zmsk == 0. ) THEN
495                     zcoff = 1.
496                  ELSE
497                     zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   &
498                        / (zmsk * 100.)
499                  ENDIF
500                  ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj)
501               ENDIF
502            END DO
503         END DO
504      ENDIF
505     
506      ! Lateral boundary conditions on ( ahm1, ahm2 )
507      !                                ==============
508      CALL lbc_lnk( ahm1, 'T', 1. )   ! T-point, unchanged sign
509      CALL lbc_lnk( ahm2, 'F', 1. )   ! F-point, unchanged sign
510
511      ! Control print
512      IF( lwp .AND. ld_print ) THEN
513         WRITE(numout,*)
514         WRITE(numout,*) 'inildf: 2D ahm1 array'
515         CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
516         WRITE(numout,*)
517         WRITE(numout,*) 'inildf: 2D ahm2 array'
518         CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
519      ENDIF
520      !
521      CALL wrk_dealloc( jpi   , jpj   , icof  )
522      !
523   END SUBROUTINE ldf_dyn_c2d_orca_R1
Note: See TracBrowser for help on using the repository browser.