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 @ 494

Last change on this file since 494 was 494, checked in by opalod, 18 years ago

nemo_v1_update_062:CE+RB: use IOM for passive tracers

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 34.1 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   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 "passivetrc_substitute.h90"
37   !!----------------------------------------------------------------------
38   !!   TOP 1.0 , LOCEAN-IPSL (2005)
39   !! $Header$
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) = 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) = 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) = 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.