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_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/TRP – NEMO

source: branches/DEV_r1879_FCM/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90 @ 2013

Last change on this file since 2013 was 2013, checked in by smasson, 14 years ago

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