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

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

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