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.
trcdmp.F90 in trunk/NEMO/TOP_SRC/TRP – NEMO

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

Last change on this file since 202 was 202, checked in by opalod, 19 years ago

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