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

Last change on this file since 367 was 352, checked in by opalod, 18 years ago

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