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 branches/dev_001_GM/NEMO/TOP_SRC/TRP – NEMO

source: branches/dev_001_GM/NEMO/TOP_SRC/TRP/trcdmp.F90 @ 772

Last change on this file since 772 was 772, checked in by gm, 16 years ago

dev_001_GM - change the name of cpp key to key_top, key_lobster, key_pisces, key_kriest and the corresponding lk_

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