source: trunk/NEMO/TOP_SRC/TRP/trcdmp.F90 @ 941

Last change on this file since 941 was 941, checked in by cetlod, 13 years ago

phasing the passive tracer transport module to the new version of NEMO, see ticket 143

  • Property svn:executable set to *
File size: 34.2 KB
Line 
1MODULE trcdmp
2   !!======================================================================
3   !!                       ***  MODULE  trcdmp  ***
4   !! Ocean physics: internal restoring trend on passive tracers
5   !!======================================================================
6#if  defined key_top  &&  defined key_trcdmp 
7   !!----------------------------------------------------------------------
8   !!   'key_top'                                                TOP models
9   !!   'key_trcdmp'                                       internal damping
10   !!----------------------------------------------------------------------
11   !!   trc_dmp      : update the tracer trend with the internal damping
12   !!   trc_dmp_init : initialization, namlist read, parameters control
13   !!   trccof_zoom  : restoring coefficient for zoom domain
14   !!   trccof       : restoring coefficient for global domain
15   !!   cofdis       : compute the distance to the coastline
16   !!----------------------------------------------------------------------
17   USE oce_trc         ! ocean dynamics and tracers variables
18   USE trc             ! ocean passive tracers variables
19   USE trctrp_lec      ! passive tracers transport
20   USE trcdta
21   USE prtctl_trc      ! Print control for debbuging
22
23   IMPLICIT NONE
24   PRIVATE
25
26   !! * Routine accessibility
27   PUBLIC trc_dmp   ! routine called by step.F90
28
29   !! * Shared module variables
30   LOGICAL , PUBLIC, PARAMETER ::   lk_trcdmp = .TRUE.    !: internal damping flag
31
32   REAL(wp), DIMENSION(jpi,jpj,jpk,jptra) ::   &
33      restotr         ! restoring coeff. on tracers (s-1)
34
35   !! * Substitutions
36#  include "top_substitute.h90"
37   !!----------------------------------------------------------------------
38   !!   TOP 1.0 , LOCEAN-IPSL (2005)
39   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/TOP_SRC/TRP/trcdmp.F90,v 1.12 2007/10/12 09:26:30 opalod Exp $
40   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
41   !!----------------------------------------------------------------------
42
43CONTAINS
44
45   SUBROUTINE trc_dmp( kt )
46      !!----------------------------------------------------------------------
47      !!                   ***  ROUTINE trc_dmp  ***
48      !!                 
49      !! ** Purpose :   Compute the passive tracer trend due to a newtonian damping
50      !!      of the tracer field towards given data field and add it to the
51      !!      general tracer trends.
52      !!
53      !! ** Method  :   Newtonian damping towards trdta computed
54      !!      and add to the general tracer trends:
55      !!                     trn = tra + restotr * (trdta - trb)
56      !!         The trend is computed either throughout the water column
57      !!      (nlmdmptr=0) or in area of weak vertical mixing (nlmdmptr=1) or
58      !!      below the well mixed layer (nlmdmptr=2)
59      !!
60      !! ** Action  : - update the tracer trends tra with the newtonian
61      !!                damping trends.
62      !!              - save the trends in trtrd ('key_trc_diatrd')
63      !!
64      !! History :
65      !!   7.0  !         (G. Madec)  Original code
66      !!        !  96-01  (G. Madec)
67      !!        !  97-05  (H. Loukos)  adapted for passive tracers
68      !!   8.5  !  02-08  (G. Madec )  free form + modules
69      !!   9.0  !  04-03  (C. Ethe)    free form + modules
70      !!----------------------------------------------------------------------
71      !! * Arguments
72      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
73
74      !! * Local declarations
75      INTEGER  ::   ji, jj, jk, jn     ! dummy loop indices
76      REAL(wp) ::   ztest, ztra, zdt   ! temporary scalars
77      CHARACTER (len=22) :: charout
78      !!----------------------------------------------------------------------
79
80      ! 0. Initialization (first time-step only)
81      !    --------------
82      IF( kt == nittrc000 ) CALL trc_dmp_init
83
84      ! 1. Newtonian damping trends on tracer fields
85      ! --------------------------------------------
86      !    compute the newtonian damping trends depending on nmldmptr
87
88!!!      zdt  = rdt * FLOAT( ndttrc )
89
90      ! Initialize the input fields for newtonian damping
91      CALL dta_trc( kt )
92
93      DO jn = 1, jptra
94
95         IF( lutini(jn) ) THEN
96
97            SELECT CASE ( nmldmptr )
98
99            CASE( 0 )                ! newtonian damping throughout the water column
100
101               DO jk = 1, jpkm1
102                  DO jj = 2, jpjm1
103                     DO ji = fs_2, fs_jpim1   ! vector opt.
104                        ztra = restotr(ji,jj,jk,jn) * ( trdta(ji,jj,jk,jn) - trb(ji,jj,jk,jn) )
105                        ! add the trends to the general tracer trends
106!!                        trn(ji,jj,jk,jn) = trn(ji,jj,jk,jn) + ztra * zdt
107                        tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
108#    if defined key_trc_diatrd
109                        ! save the trends for diagnostics
110                        IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),jpdiatrc-1) = ztra
111#    endif
112                     END DO
113                  END DO
114               END DO
115
116            CASE ( 1 )                ! no damping in the turbocline (avt > 5 cm2/s)
117               DO jk = 1, jpkm1
118                  DO jj = 2, jpjm1
119                     DO ji = fs_2, fs_jpim1   ! vector opt.
120                        ztest = avt(ji,jj,jk) - 5.e-4
121                        IF( ztest < 0. ) THEN
122                           ztra = restotr(ji,jj,jk,jn) * ( trdta(ji,jj,jk,jn) - trb(ji,jj,jk,jn) )
123                        ELSE
124                           ztra = 0.e0
125                        ENDIF
126                        ! add the trends to the general tracer trends
127!!                        trn(ji,jj,jk,jn) = trn(ji,jj,jk,jn) + ztra * zdt
128                        tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra 
129#    if defined key_trc_diatrd
130                        ! save the trends for diagnostics
131                        IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),jpdiatrc-1) = ztra
132#    endif
133                     END DO
134                  END DO
135               END DO
136
137            CASE ( 2 )                ! no damping in the mixed layer
138               DO jk = 1, jpkm1
139                  DO jj = 2, jpjm1
140                     DO ji = fs_2, fs_jpim1   ! vector opt.
141                        IF( fsdept(ji,jj,jk) >= hmlp (ji,jj) ) THEN
142                           ztra = restotr(ji,jj,jk,jn) * ( trdta(ji,jj,jk,jn) - trb(ji,jj,jk,jn) )
143                        ELSE
144                           ztra = 0.e0
145                        ENDIF
146                        ! add the trends to the general tracer trends
147!!                        trn(ji,jj,jk,jn) = trn(ji,jj,jk,jn) + ztra * zdt
148                        tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
149#    if defined key_trc_diatrd
150                        ! save the trends for diagnostics
151                        IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),jpdiatrc-1) = ztra
152#    endif
153                     END DO
154                  END DO
155               END DO
156               
157            END SELECT
158
159         ENDIF
160
161      END DO
162
163     IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
164         WRITE(charout, FMT="('dmp')")
165         CALL prt_ctl_trc_info(charout)
166         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
167      ENDIF
168 
169      trb(:,:,:,:) = trn(:,:,:,:)
170   
171   END SUBROUTINE trc_dmp
172
173
174   SUBROUTINE trc_dmp_init
175      !!----------------------------------------------------------------------
176      !!                  ***  ROUTINE trc_dmp_init  ***
177      !!
178      !! ** Purpose :   Initialization for the newtonian damping
179      !!
180      !! ** Method  :   read the nammbf namelist and check the parameters
181      !!      called by trc_dmp at the first timestep (nit000)
182      !!
183      !! History :
184      !!   8.5  !  02-08  (G. Madec)  Original code
185      !!----------------------------------------------------------------------
186
187      SELECT CASE ( ndmptr )
188
189      CASE ( -1 )               ! ORCA: damping in Red & Med Seas only
190         IF(lwp) WRITE(numout,*) '          tracer damping in the Med & Red seas only'
191
192      CASE ( 1:90 )             ! Damping poleward of 'ndmptr' degrees
193         IF(lwp) WRITE(numout,*) '          tracer damping poleward of', ndmptr, ' degrees'
194
195      CASE DEFAULT
196         WRITE(ctmp1,*) '          bad flag value for ndmp = ', ndmp
197         CALL ctl_stop(ctmp1)
198
199      END SELECT
200
201
202      SELECT CASE ( nmldmptr )
203
204      CASE ( 0 )                ! newtonian damping throughout the water column
205         IF(lwp) WRITE(numout,*) '          tracer damping throughout the water column'
206
207      CASE ( 1 )                ! no damping in the turbocline (avt > 5 cm2/s)
208         IF(lwp) WRITE(numout,*) '          no tracer damping in the turbocline'
209
210      CASE ( 2 )                ! no damping in the mixed layer
211         IF(lwp) WRITE(numout,*) '          no tracer damping in the mixed layer'
212
213      CASE DEFAULT
214         WRITE(ctmp1,*) '          bad flag value for nmldmp = ', nmldmp
215         CALL ctl_stop(ctmp1)
216
217
218      END SELECT
219
220
221      ! 3. Damping coefficients initialization
222     ! --------------------------------------
223
224         IF( lzoom ) THEN
225            CALL trccof_zoom
226         ELSE
227            CALL trccof
228         ENDIF
229 
230   END SUBROUTINE trc_dmp_init
231
232
233   SUBROUTINE trccof_zoom
234      !!----------------------------------------------------------------------
235      !!                  ***  ROUTINE trccof_zoom  ***
236      !!
237      !! ** Purpose :   Compute the damping coefficient for zoom domain
238      !!
239      !! ** Method  : - set along closed boundary due to zoom a damping over
240      !!      6 points with a max time scale of 5 days.
241      !!              - ORCA arctic/antarctic zoom: set the damping along
242      !!      south/north boundary over a latitude strip.
243      !!
244      !! ** Action  : - restotr, the damping coeff. passive tracers
245      !!
246      !! History :
247      !!   9.0  !  03-09  (G. Madec)  Original code
248      !!   9.0  !  04-03  (C. Ethe)   adapted for passive tracers
249      !!----------------------------------------------------------------------
250      !! * Local declarations
251      INTEGER ::   ji, jj, jk, jn       ! dummy loop indices
252      REAL(wp) ::   &
253         zlat, zlat0, zlat1, zlat2     ! temporary scalar
254      REAL(wp), DIMENSION(6)  ::   &
255         zfact                         ! temporary workspace
256      !!----------------------------------------------------------------------
257
258      zfact(1) =  1.
259      zfact(2) =  1. 
260      zfact(3) = 11./12.
261      zfact(4) =  8./12.
262      zfact(5) =  4./12.
263      zfact(6) =  1./12.
264      zfact(:) = zfact(:) / ( 5. * rday )    ! 5 days max restoring time scale
265
266      restotr(:,:,:,:) = 0.e0
267
268      ! damping along the forced closed boundary over 6 grid-points
269      DO jn = 1, 6
270         IF( lzoom_w )   restotr( mi0(jn+jpizoom):mi1(jn+jpizoom), : , : , : ) = zfact(jn) ! west  closed
271         IF( lzoom_s )   restotr( : , mj0(jn+jpjzoom):mj1(jn+jpjzoom), : , : ) = zfact(jn) ! south closed
272         IF( lzoom_e )   restotr( mi0(jpiglo+jpizoom-1-jn):mi1(jpiglo+jpizoom-1-jn) , : , : , : ) &
273                       &              = zfact(jn)                                 ! east  closed
274         IF( lzoom_n )   restotr( : , mj0(jpjglo+jpjzoom-1-jn):mj1(jpjglo+jpjzoom-1-jn) , : , : ) &
275                       &              = zfact(jn)                                 ! north closed
276      END DO
277
278
279      IF( lzoom_arct .AND. lzoom_anta ) THEN
280
281         ! ====================================================
282         !  ORCA configuration : arctic zoom or antarctic zoom
283         ! ====================================================
284
285         IF(lwp) WRITE(numout,*)
286         IF(lwp .AND. lzoom_arct ) WRITE(numout,*) '              trccof_zoom : ORCA    Arctic zoom'
287         IF(lwp .AND. lzoom_arct ) WRITE(numout,*) '              trccof_zoom : ORCA Antarctic zoom'
288         IF(lwp) WRITE(numout,*)
289
290         ! ... Initialization :
291         !     zlat0 : latitude strip where resto decreases
292         !     zlat1 : resto = 1 before zlat1
293         !     zlat2 : resto decreases from 1 to 0 between zlat1 and zlat2
294         restotr(:,:,:,:) = 0.e0
295         zlat0 = 10.
296         zlat1 = 30.
297         zlat2 = zlat1 + zlat0
298
299         ! ... Compute arrays resto ; value for internal damping : 5 days
300         DO jn = 1, jptra
301            DO jk = 2, jpkm1 
302               DO jj = 1, jpj
303                  DO ji = 1, jpi
304                     zlat = ABS( gphit(ji,jj) )
305                     IF ( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN
306                        restotr(ji,jj,jk,jn) = 0.5 * ( 1./(5.*rday) ) *   &
307                           ( 1. - COS(rpi*(zlat2-zlat)/zlat0) ) 
308                     ELSE IF ( zlat < zlat1 ) THEN
309                        restotr(ji,jj,jk,jn) = 1./(5.*rday)
310                     ENDIF
311                  END DO
312               END DO
313            END DO
314         END DO
315
316      ENDIF
317
318      ! ... Mask resto array
319        DO jn = 1, jptra
320           restotr(:,:,:,jn) = restotr(:,:,:,jn) * tmask(:,:,:)
321        END DO
322
323
324   END SUBROUTINE trccof_zoom
325
326   SUBROUTINE trccof
327      !!----------------------------------------------------------------------
328      !!                  ***  ROUTINE trccof  ***
329      !!
330      !! ** Purpose :   Compute the damping coefficient
331      !!
332      !! ** Method  :   Arrays defining the damping are computed for each grid
333      !!      point passive tracers (restotr)
334      !!      Damping depends on distance to coast, depth and latitude
335      !!
336      !! ** Action  : - restotr, the damping coeff. for passive tracers
337      !!
338      !! History :
339      !!   5.0  !  91-03  (O. Marti, G. Madec)  Original code
340      !!        !  92-06  (M. Imbard)  doctor norme
341      !!        !  96-01  (G. Madec) statement function for e3
342      !!        !  98-07  (M. Imbard, G. Madec) ORCA version
343      !!        !  00-08  (G. Madec, D. Ludicone)
344      !!   8.2  !  04-03  (H. Loukos) adapted for passive tracers
345      !!        !  04-02  (O. Aumont, C. Ethe) rewritten for debuging and update
346      !!----------------------------------------------------------------------
347      !! * Modules used
348      USE iom
349      USE ioipsl
350
351      !! * Local declarations
352      INTEGER ::  ji, jj, jk, jn    ! dummy loop indices
353      INTEGER ::   itime
354      INTEGER ::  ii0, ii1, ij0, ij1  !    "          "
355      INTEGER ::   &
356         idmp,     &  ! logical unit for file restoring damping term
357         icot         ! logical unit for file distance to the coast
358
359      CHARACTER (len=32) ::  clname, clname2, clname3
360      REAL(wp) ::   &
361         zdate0, zinfl, zlon,         & ! temporary scalars
362         zlat, zlat0, zlat1, zlat2,   & !    "         "
363         zsdmp, zbdmp                   !    "         "
364      REAL(wp), DIMENSION(jpk) ::   &
365         gdept, zhfac
366      REAL(wp), DIMENSION(jpi,jpj) ::   &
367         zmrs
368      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
369         zdct
370      !!----------------------------------------------------------------------
371
372      ! ====================================
373      !  ORCA configuration : global domain
374      ! ====================================
375
376      IF(lwp) WRITE(numout,*)
377      IF(lwp) WRITE(numout,*) '              trccof : Global domain of ORCA'
378      IF(lwp) WRITE(numout,*) '              ------------------------------'
379
380
381      ! ... Initialization :
382      !   zdct()      : distant to the coastline
383      !   resto()     : array of restoring coeff.
384     
385      zdct (:,:,:) = 0.e0
386      restotr(:,:,:,:) = 0.e0
387
388
389      IF ( ndmptr > 0 ) THEN
390
391         !    ------------------------------------
392         !     Damping poleward of 'ndmptr' degrees
393         !    ------------------------------------
394
395         IF(lwp) WRITE(numout,*)
396         IF(lwp) WRITE(numout,*) '              Damping poleward of ', ndmptr,' deg.'
397         IF(lwp) WRITE(numout,*)
398
399         ! ... Distance to coast (zdct)
400
401         IF(lwp) WRITE(numout,*)
402         IF(lwp) WRITE(numout,*) ' dtacof : distance to coast file'
403         CALL iom_open ( 'dist.coast.trc.nc', icot )
404         IF( icot > 0 ) THEN
405            CALL iom_get ( icot, jpdom_data, 'Tcoast', zdct )
406            CALL iom_close (icot)
407         ELSE
408            !   ... Compute and save the distance-to-coast array (output in zdct)
409            CALL cofdis( zdct )
410         ENDIF
411
412
413         ! ... Compute arrays resto
414         !      zinfl : distance of influence for damping term
415         !      zlat0 : latitude strip where resto decreases
416         !      zlat1 : resto = 0 between -zlat1 and zlat1
417         !      zlat2 : resto increases from 0 to 1 between |zlat1| and |zlat2|
418         !          and resto = 1 between |zlat2| and 90 deg.
419         zinfl = 1000.e3
420         zlat0 = 10
421         zlat1 = ndmptr
422         zlat2 = zlat1 + zlat0
423
424         DO jn = 1, jptra
425            DO jj = 1, jpj
426               DO ji = 1, jpi
427                  zlat = ABS( gphit(ji,jj) )
428                  IF ( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN
429                     restotr(ji,jj,1,jn) = 0.5 * ( 1. - COS(rpi*(zlat-zlat1)/zlat0 ) )
430                  ELSEIF ( zlat > zlat2 ) THEN
431                     restotr(ji,jj,1,jn) = 1.
432                  ENDIF
433               END DO
434            END DO
435         END DO
436
437         !   ... North Indian ocean (20N/30N x 45E/100E) : resto=0
438         IF ( ndmptr == 20 ) THEN
439            DO jn = 1, jptra
440               DO jj = 1, jpj
441                  DO ji = 1, jpi
442                     zlat = gphit(ji,jj)
443                     zlon = MOD( glamt(ji,jj), 360. )
444                     IF ( zlat1 < zlat .AND. zlat < zlat2 .AND.   &
445                        45.  < zlon .AND. zlon < 100. ) THEN
446                        restotr(ji,jj,1,jn) = 0.
447                     ENDIF
448                  END DO
449               END DO
450            END DO
451         ENDIF
452
453         zsdmp = 1./(sdmptr * rday)
454         zbdmp = 1./(bdmptr * rday)
455         DO jn = 1, jptra
456            DO jk = 2, jpkm1
457               DO jj = 1, jpj
458                  DO ji = 1, jpi
459                     zdct(ji,jj,jk) = MIN( zinfl, zdct(ji,jj,jk) )
460
461                     !   ... Decrease the value in the vicinity of the coast
462                     restotr(ji,jj,jk,jn) = restotr(ji,jj,1,jn)*0.5   &
463                        &                 * ( 1. - COS( rpi*zdct(ji,jj,jk)/zinfl) )
464
465                     !   ... Vertical variation from zsdmp (sea surface) to zbdmp (bottom)
466                     restotr(ji,jj,jk,jn) = restotr(ji,jj,jk,jn)   &
467                        &                 * ( zbdmp + (zsdmp-zbdmp)*EXP(-fsdept(ji,jj,jk)/hdmptr) )
468                  END DO
469               END DO
470            END DO
471         END DO
472
473      ENDIF
474
475
476      IF( cp_cfg == "orca" .AND. ( ndmptr > 0 .OR. ndmptr == -1 ) ) THEN
477
478         !                                         ! =========================
479         !                                         !  Med and Red Sea damping
480         !                                         ! =========================
481         IF(lwp)WRITE(numout,*)
482         IF(lwp)WRITE(numout,*) '              ORCA configuration: Damping in Med and Red Seas'
483
484
485         zmrs(:,:) = 0.e0                             ! damping term on the Med or Red Sea
486
487         SELECT CASE ( jp_cfg )
488            !                                           ! =======================
489         CASE ( 4 )                                     !  ORCA_R4 configuration
490            !                                           ! =======================
491
492            ! Mediterranean Sea
493            ij0 =  50   ;   ij1 =  56
494            ii0 =  81   ;   ii1 =  91   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
495            ij0 =  50   ;   ij1 =  55
496            ii0 =  75   ;   ii1 =  80   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
497            ij0 =  52   ;   ij1 =  53
498            ii0 =  70   ;   ii1 =  74   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
499            ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea
500            DO jk = 1, 17
501               zhfac (jk) = 0.5*( 1.- COS( rpi*(jk-1)/16. ) ) / rday
502            END DO
503            DO jk = 18, jpkm1
504               zhfac (jk) = 1./rday
505            END DO
506
507            !                                        ! =======================
508         CASE ( 2 )                                  !  ORCA_R2 configuration
509            !                                        ! =======================
510
511            ! Mediterranean Sea
512            ij0 =  96   ;   ij1 = 110
513            ii0 = 157   ;   ii1 = 181   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
514            ij0 = 100   ;   ij1 = 110
515            ii0 = 144   ;   ii1 = 156   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
516            ij0 = 100   ;   ij1 = 103
517            ii0 = 139   ;   ii1 = 143   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
518            ! Decrease before Gibraltar Strait
519            ij0 = 101   ;   ij1 = 102
520            ii0 = 139   ;   ii1 = 141   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.e0
521            ii0 = 142   ;   ii1 = 142   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0 / 90.e0
522            ii0 = 143   ;   ii1 = 143   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40e0
523            ii0 = 144   ;   ii1 = 144   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75e0
524            ! Red Sea
525            ij0 =  87   ;   ij1 =  96
526            ii0 = 147   ;   ii1 = 163   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
527            ! Decrease before Bab el Mandeb Strait
528            ij0 =  91   ;   ij1 =  91
529            ii0 = 153   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.80e0
530            ij0 =  90   ;   ij1 =  90
531            ii0 = 153   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40e0
532            ij0 =  89   ;   ij1 =  89
533            ii0 = 158   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0 / 90.e0
534            ij0 =  88   ;   ij1 =  88
535            ii0 = 160   ;   ii1 = 163   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.e0
536            ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea
537            DO jk = 1, 17
538               zhfac (jk) = 0.5*( 1.- COS( rpi*(jk-1)/16. ) ) / rday
539            END DO
540            DO jk = 18, jpkm1
541               zhfac (jk) = 1./rday
542            END DO
543
544            !                                        ! =======================
545         CASE ( 05 )                                 !  ORCA_R05 configuration
546            !                                        ! =======================
547
548            ! Mediterranean Sea
549            ii0 = 568   ;   ii1 = 574 
550            ij0 = 324   ;   ij1 = 333   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
551            ii0 = 575   ;   ii1 = 658
552            ij0 = 314   ;   ij1 = 366   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
553            ! Black Sea (remaining part
554            ii0 = 641   ;   ii1 = 651
555            ij0 = 367   ;   ij1 = 372   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
556            ! Decrease before Gibraltar Strait
557            ii0 = 324   ;   ii1 = 333
558            ij0 = 565   ;   ij1 = 565   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0 / 90.e0
559            ij0 = 566   ;   ij1 = 566   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
560            ij0 = 567   ;   ij1 = 567   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75
561            ! Red Sea
562            ii0 = 641   ;   ii1 = 665
563            ij0 = 270   ;   ij1 = 310   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
564            ! Decrease before Bab el Mandeb Strait
565            ii0 = 666   ;   ii1 = 675
566            ij0 = 270   ;   ij1 = 290   
567            DO ji = mi0(ii0), mi1(ii1)
568               zmrs( ji , mj0(ij0):mj1(ij1) ) = 0.1 * ABS( FLOAT(ji - mi1(ii1)) )
569            END DO
570            zsdmp = 1./(sdmptr * rday)
571            zbdmp = 1./(bdmptr * rday)
572            DO jk = 1, jpk
573               zhfac (jk) = ( zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(1,1,jk)/hdmptr) )
574            END DO
575
576            !                                       ! ========================
577         CASE ( 025 )                               !  ORCA_R025 configuration
578
579            CALL ctl_stop( ' Not yet implemented in ORCA_R025' ) 
580
581         END SELECT
582
583         DO jn = 1, jptra
584            DO jk = 1, jpkm1
585               restotr(:,:,jk,jn) = zmrs(:,:) * zhfac(jk) + ( 1. - zmrs(:,:) ) * restotr(:,:,jk,jn)
586            END DO
587
588            ! Mask resto array and set to 0 first and last levels
589            restotr(:,:, : ,jn) = restotr(:,:,:,jn) * tmask(:,:,:)
590            restotr(:,:, 1 ,jn) = 0.e0
591            restotr(:,:,jpk,jn) = 0.e0
592         END DO
593
594      ELSE
595         !    ------------
596         !     No damping
597         !    ------------
598         CALL ctl_stop( 'Choose a correct value of ndmp or DO NOT defined key_tradmp' )
599
600      ENDIF
601
602        !    ----------------------------
603         !     Create Print damping array
604         !    ----------------------------
605         
606         ! ndmpftr   : = 1 create a damping.coeff NetCDF file
607
608      IF( ndmpftr == 1 ) THEN
609         DO jn = 1, jptra
610            IF(lwp) WRITE(numout,*) '  create damping.coeff.nc file  ',jn
611            itime   = 0
612            clname3 = 'damping.coeff'//ctrcnm(jn)
613            CALL ymds2ju( 0     , 1     , 1      , 0.e0 , zdate0 )
614            CALL restini( 'NONE', jpi   , jpj    , glamt, gphit,    &
615           &              jpk   , gdept , clname3, itime, zdate0,   &
616           &              rdt   , idmp  , domain_id=nidom)
617            CALL restput( idmp, 'Resto', jpi, jpj, jpk, 0 , restotr(:,:,:,jn)  )
618            CALL restclo( idmp )
619         END DO
620      ENDIF
621
622
623   END SUBROUTINE trccof
624
625
626   SUBROUTINE cofdis ( pdct )
627      !!----------------------------------------------------------------------
628      !!                 ***  ROUTINE cofdis  ***
629      !!
630      !! ** Purpose :   Compute the distance between ocean T-points and the
631      !!      ocean model coastlines. Save the distance in a NetCDF file.
632      !!
633      !! ** Method  :   For each model level, the distance-to-coast is
634      !!      computed as follows :
635      !!       - The coastline is defined as the serie of U-,V-,F-points
636      !!      that are at the ocean-land bound.
637      !!       - For each ocean T-point, the distance-to-coast is then
638      !!      computed as the smallest distance (on the sphere) between the
639      !!      T-point and all the coastline points.
640      !!       - For land T-points, the distance-to-coast is set to zero.
641      !!      C A U T I O N : Computation not yet implemented in mpp case.
642      !!
643      !! ** Action  : - pdct, distance to the coastline (argument)
644      !!              - NetCDF file 'trc.dist.coast.nc'
645      !!       
646      !! History :
647      !!   7.0  !  01-02  (M. Imbard)  Original code
648      !!   8.1  !  01-02  (G. Madec, E. Durand)
649      !!   8.5  !  02-08  (G. Madec, E. Durand)  Free form, F90
650      !!----------------------------------------------------------------------
651      !! * Modules used
652      USE ioipsl
653
654      !! * Arguments
655      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) ::   &
656         pdct                     ! distance to the coastline
657
658      !! * local declarations
659      INTEGER :: ji, jj, jk, jl      ! dummy loop indices
660      INTEGER :: iju, ijt            ! temporary integers
661      INTEGER :: icoast, itime
662      INTEGER ::   &
663         icot         ! logical unit for file distance to the coast
664      LOGICAL, DIMENSION(jpi,jpj) ::   &
665         llcotu, llcotv, llcotf   ! ???
666      CHARACTER (len=32) ::   clname
667      REAL(wp) ::   zdate0
668      REAL(wp), DIMENSION(jpi,jpj) ::   &
669         zxt, zyt, zzt,                 &  ! cartesian coordinates for T-points
670         zmask                             
671      REAL(wp), DIMENSION(3*jpi*jpj) ::   &
672         zxc, zyc, zzc, zdis      ! temporary workspace
673      !!----------------------------------------------------------------------
674
675      ! 0. Initialization
676      ! -----------------
677      IF(lwp) WRITE(numout,*)
678      IF(lwp) WRITE(numout,*) 'cofdis : compute the distance to coastline'
679      IF(lwp) WRITE(numout,*) '~~~~~~'
680      IF(lwp) WRITE(numout,*)
681      IF( lk_mpp ) &
682           & CALL ctl_stop('         Computation not yet implemented with key_mpp_...', &
683           &               '         Rerun the code on another computer or ', &
684           &               '         create the "dist.coast.nc" file using IDL' )
685
686
687      pdct(:,:,:) = 0.e0
688      zxt(:,:) = cos( rad * gphit(:,:) ) * cos( rad * glamt(:,:) )
689      zyt(:,:) = cos( rad * gphit(:,:) ) * sin( rad * glamt(:,:) )
690      zzt(:,:) = sin( rad * gphit(:,:) )
691
692
693      ! 1. Loop on vertical levels
694      ! --------------------------
695      !                                                ! ===============
696      DO jk = 1, jpkm1                                 ! Horizontal slab
697         !                                             ! ===============
698         ! Define the coastline points (U, V and F)
699         DO jj = 2, jpjm1
700            DO ji = 2, jpim1
701               zmask(ji,jj) =  ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
702                   &           + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) )
703               llcotu(ji,jj) = ( tmask(ji,jj,  jk) + tmask(ji+1,jj  ,jk) == 1. ) 
704               llcotv(ji,jj) = ( tmask(ji,jj  ,jk) + tmask(ji  ,jj+1,jk) == 1. ) 
705               llcotf(ji,jj) = ( zmask(ji,jj) > 0. ) .AND. ( zmask(ji,jj) < 4. )
706            END DO
707         END DO
708
709         ! Lateral boundaries conditions
710         llcotu(:, 1 ) = umask(:,  2  ,jk) == 1
711         llcotu(:,jpj) = umask(:,jpjm1,jk) == 1
712         llcotv(:, 1 ) = vmask(:,  2  ,jk) == 1
713         llcotv(:,jpj) = vmask(:,jpjm1,jk) == 1
714         llcotf(:, 1 ) = fmask(:,  2  ,jk) == 1
715         llcotf(:,jpj) = fmask(:,jpjm1,jk) == 1
716
717         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
718            llcotu( 1 ,:) = llcotu(jpim1,:)
719            llcotu(jpi,:) = llcotu(  2  ,:)
720            llcotv( 1 ,:) = llcotv(jpim1,:)
721            llcotv(jpi,:) = llcotv(  2  ,:)
722            llcotf( 1 ,:) = llcotf(jpim1,:)
723            llcotf(jpi,:) = llcotf(  2  ,:)
724         ELSE
725            llcotu( 1 ,:) = umask(  2  ,:,jk) == 1
726            llcotu(jpi,:) = umask(jpim1,:,jk) == 1
727            llcotv( 1 ,:) = vmask(  2  ,:,jk) == 1
728            llcotv(jpi,:) = vmask(jpim1,:,jk) == 1
729            llcotf( 1 ,:) = fmask(  2  ,:,jk) == 1
730            llcotf(jpi,:) = fmask(jpim1,:,jk) == 1
731         ENDIF
732         IF( nperio == 3 .OR. nperio == 4 ) THEN
733            DO ji = 1, jpim1
734               iju = jpi - ji + 1
735               llcotu(ji,jpj  ) = llcotu(iju,jpj-2)
736               llcotf(ji,jpj-1) = llcotf(iju,jpj-2)
737               llcotf(ji,jpj  ) = llcotf(iju,jpj-3)
738            END DO
739            DO ji = jpi/2, jpi-1
740               iju = jpi - ji + 1
741               llcotu(ji,jpjm1) = llcotu(iju,jpjm1)
742            END DO
743            DO ji = 2, jpi
744               ijt = jpi - ji + 2
745               llcotv(ji,jpj-1) = llcotv(ijt,jpj-2)
746               llcotv(ji,jpj  ) = llcotv(ijt,jpj-3)
747            END DO
748         ENDIF
749         IF( nperio == 5 .OR. nperio == 6 ) THEN
750            DO ji = 1, jpim1
751               iju = jpi - ji
752               llcotu(ji,jpj  ) = llcotu(iju,jpj-1)
753               llcotf(ji,jpj  ) = llcotf(iju,jpj-2)
754            END DO
755            DO ji = jpi/2, jpi-1
756               iju = jpi - ji
757               llcotf(ji,jpjm1) = llcotf(iju,jpjm1)
758            END DO
759            DO ji = 1, jpi
760               ijt = jpi - ji + 1
761               llcotv(ji,jpj  ) = llcotv(ijt,jpj-1)
762            END DO
763            DO ji = jpi/2+1, jpi
764               ijt = jpi - ji + 1
765               llcotv(ji,jpjm1) = llcotv(ijt,jpjm1)
766            END DO
767         ENDIF
768
769         ! Compute cartesian coordinates of coastline points
770         ! and the number of coastline points
771
772         icoast = 0
773         DO jj = 1, jpj
774            DO ji = 1, jpi
775               IF( llcotf(ji,jj) ) THEN
776                  icoast = icoast + 1
777                  zxc(icoast) = COS( rad*gphif(ji,jj) ) * COS( rad*glamf(ji,jj) )
778                  zyc(icoast) = COS( rad*gphif(ji,jj) ) * SIN( rad*glamf(ji,jj) )
779                  zzc(icoast) = SIN( rad*gphif(ji,jj) )
780               ENDIF
781               IF( llcotu(ji,jj) ) THEN
782                  icoast = icoast+1
783                  zxc(icoast) = COS( rad*gphiu(ji,jj) ) * COS( rad*glamu(ji,jj) )
784                  zyc(icoast) = COS( rad*gphiu(ji,jj) ) * SIN( rad*glamu(ji,jj) )
785                  zzc(icoast) = SIN( rad*gphiu(ji,jj) )
786               ENDIF
787               IF( llcotv(ji,jj) ) THEN
788                  icoast = icoast+1
789                  zxc(icoast) = COS( rad*gphiv(ji,jj) ) * COS( rad*glamv(ji,jj) )
790                  zyc(icoast) = COS( rad*gphiv(ji,jj) ) * SIN( rad*glamv(ji,jj) )
791                  zzc(icoast) = SIN( rad*gphiv(ji,jj) )
792               ENDIF
793            END DO
794         END DO
795
796         ! Distance for the T-points
797
798         DO jj = 1, jpj
799            DO ji = 1, jpi
800               IF( tmask(ji,jj,jk) == 0. ) THEN
801                  pdct(ji,jj,jk) = 0.
802               ELSE
803                  DO jl = 1, icoast
804                     zdis(jl) = ( zxt(ji,jj) - zxc(jl) )**2   &
805                              + ( zyt(ji,jj) - zyc(jl) )**2   &
806                              + ( zzt(ji,jj) - zzc(jl) )**2
807                  END DO
808                  pdct(ji,jj,jk) = ra * SQRT( MINVAL( zdis(1:icoast) ) )
809               ENDIF
810            END DO
811         END DO
812         !                                                ! ===============
813      END DO                                              !   End of slab
814      !                                                   ! ===============
815
816
817      ! 2. Create the  distance to the coast file in NetCDF format
818      ! ----------------------------------------------------------   
819      clname = 'trc.dist.coast'
820      itime = 0
821      CALL ymds2ju( 0     , 1     , 1     , 0.e0 , zdate0 )
822      CALL restini( 'NONE', jpi   , jpj   , glamt, gphit ,   &
823                    jpk   , gdept , clname, itime, zdate0,   &
824                    rdt   , icot , domain_id=nidom         )
825      CALL restput( icot, 'Tcoast', jpi, jpj, jpk, 0, pdct )
826      CALL restclo( icot )
827
828   END SUBROUTINE cofdis
829
830#else
831   !!----------------------------------------------------------------------
832   !!   Default key                                     NO internal damping
833   !!----------------------------------------------------------------------
834   LOGICAL , PUBLIC, PARAMETER ::   lk_trcdmp = .FALSE.    !: internal damping flag
835CONTAINS
836   SUBROUTINE trc_dmp( kt )        ! Empty routine
837      INTEGER, INTENT(in) :: kt
838      WRITE(*,*) 'trc_dmp: You should not have seen this print! error?', kt
839   END SUBROUTINE trc_dmp
840#endif
841
842   !!======================================================================
843END MODULE trcdmp
Note: See TracBrowser for help on using the repository browser.