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

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

nemo_v1_update_05:RB+OA:Update and rewritting of (part of) the TOP component

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