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.
tradmp.F90 in trunk/NEMO/OPA_SRC/TRA – NEMO

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

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

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