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

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

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