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 tags/nemo_dev_x1/NEMO/OPA_SRC/TRA – NEMO

source: tags/nemo_dev_x1/NEMO/OPA_SRC/TRA/tradmp.F90 @ 1512

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

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