source: trunk/NEMO/OPA_SRC/TRA/tradmp.F90 @ 392

Last change on this file since 392 was 392, checked in by opalod, 15 years ago

RB:nemo_v1_update_038: first integration of Agrif :

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