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 branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90 @ 2287

Last change on this file since 2287 was 2287, checked in by smasson, 14 years ago

update licence of all NEMO files...
  • Property svn:keywords set to Id
File size: 35.9 KB
Line 
1MODULE tradmp
2   !!======================================================================
3   !!                       ***  MODULE  tradmp  ***
4   !! Ocean physics: internal restoring trend on active tracers (T and S)
5   !!======================================================================
6   !! History :  OPA  ! 1991-03  (O. Marti, G. Madec)  Original code
7   !!                 ! 1992-06  (M. Imbard)  doctor norme
8   !!                 ! 1996-01  (G. Madec)  statement function for e3
9   !!                 ! 1997-05  (G. Madec)  macro-tasked on jk-slab
10   !!                 ! 1998-07  (M. Imbard, G. Madec) ORCA version
11   !!            7.0  ! 2001-02  (M. Imbard)  cofdis, Original code
12   !!            8.1  ! 2001-02  (G. Madec, E. Durand)  cleaning
13   !!  NEMO      1.0  ! 2002-08  (G. Madec, E. Durand)  free form + modules
14   !!            3.2  ! 2009-08  (G. Madec, C. Talandier)  DOCTOR norm for namelist parameter
15   !!            3.3  ! 2010-06  (C. Ethe, G. Madec) merge TRA-TRC
16   !!----------------------------------------------------------------------
17#if   defined key_tradmp   ||   defined key_esopa
18   !!----------------------------------------------------------------------
19   !!   key_tradmp                                         internal damping
20   !!----------------------------------------------------------------------
21   !!   tra_dmp      : update the tracer trend with the internal damping
22   !!   tra_dmp_init : initialization, namlist read, parameters control
23   !!   dtacof_zoom  : restoring coefficient for zoom domain
24   !!   dtacof       : restoring coefficient for global domain
25   !!   cofdis       : compute the distance to the coastline
26   !!----------------------------------------------------------------------
27   USE oce             ! ocean dynamics and tracers variables
28   USE dom_oce         ! ocean space and time domain variables
29   USE trdmod_oce         ! ocean space and time domain variables
30   USE trdtra         ! ocean space and time domain variables
31   USE zdf_oce         ! ocean vertical physics
32   USE phycst          ! Define parameters for the routines
33   USE dtatem          ! temperature data
34   USE dtasal          ! salinity data
35   USE zdfmxl          ! mixed layer depth
36   USE in_out_manager  ! I/O manager
37   USE lib_mpp         ! distribued memory computing
38   USE prtctl          ! Print control
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC   tra_dmp      ! routine called by step.F90
44   PUBLIC   tra_dmp_init ! routine called by opa.F90
45   PUBLIC   dtacof       ! routine called by tradmp.F90 and trcdmp.F90
46   PUBLIC   dtacof_zoom  ! routine called by tradmp.F90 and trcdmp.F90
47
48#if ! defined key_agrif
49   LOGICAL, PUBLIC, PARAMETER ::   lk_tradmp = .TRUE.     !: internal damping flag
50#else
51   LOGICAL, PUBLIC            ::   lk_tradmp = .TRUE.     !: internal damping flag
52#endif
53   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   strdmp   !: damping salinity trend (psu/s)
54   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   ttrdmp   !: damping temperature trend (Centigrade/s)
55   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   resto    !: restoring coeff. on T and S (s-1)
56   
57   !                             !!* Namelist namtra_dmp : T & S newtonian damping *
58   INTEGER  ::   nn_hdmp =   -1   ! = 0/-1/'latitude' for damping over T and S
59   INTEGER  ::   nn_zdmp =    0   ! = 0/1/2 flag for damping in the mixed layer
60   REAL(wp) ::   rn_surf =   50.  ! surface time scale for internal damping        [days]
61   REAL(wp) ::   rn_bot  =  360.  ! bottom time scale for internal damping         [days]
62   REAL(wp) ::   rn_dep  =  800.  ! depth of transition between rn_surf and rn_bot [meters]
63   INTEGER  ::   nn_file =    2   ! = 1 create a damping.coeff NetCDF file
64
65   !! * Substitutions
66#  include "domzgr_substitute.h90"
67#  include "vectopt_loop_substitute.h90"
68   !!----------------------------------------------------------------------
69   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
70   !! $Id$
71   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
72   !!----------------------------------------------------------------------
73
74CONTAINS
75
76   SUBROUTINE tra_dmp( kt )
77      !!----------------------------------------------------------------------
78      !!                   ***  ROUTINE tra_dmp  ***
79      !!                 
80      !! ** Purpose :   Compute the tracer trend due to a newtonian damping
81      !!      of the tracer field towards given data field and add it to the
82      !!      general tracer trends.
83      !!
84      !! ** Method  :   Newtonian damping towards t_dta and s_dta computed
85      !!      and add to the general tracer trends:
86      !!                     ta = ta + resto * (t_dta - tb)
87      !!                     sa = sa + resto * (s_dta - sb)
88      !!         The trend is computed either throughout the water column
89      !!      (nlmdmp=0) or in area of weak vertical mixing (nlmdmp=1) or
90      !!      below the well mixed layer (nlmdmp=2)
91      !!
92      !! ** Action  : - (ta,sa)   tracer trends updated with the damping trend
93      !!----------------------------------------------------------------------
94      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
95      !!
96      REAL(wp) ::   zta, zsa       ! temporary scalars
97      INTEGER ::   ji, jj, jk   ! dummy loop indices
98      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ztrdt, ztrds
99      !!----------------------------------------------------------------------
100
101      IF( l_trdtra )   THEN                    !* Save ta and sa trends
102         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
103         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsa(:,:,:,jp_sal)
104      ENDIF
105
106      SELECT CASE ( nn_zdmp ) 
107      !
108      CASE( 0 )                   !==  newtonian damping throughout the water column  ==!
109         DO jk = 1, jpkm1
110            DO jj = 2, jpjm1
111               DO ji = fs_2, fs_jpim1   ! vector opt.
112                  zta = resto(ji,jj,jk) * ( t_dta(ji,jj,jk) - tsb(ji,jj,jk,jp_tem) )
113                  zsa = resto(ji,jj,jk) * ( s_dta(ji,jj,jk) - tsb(ji,jj,jk,jp_sal) )
114                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta
115                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
116                  ! save the salinity trend (used in asmtrj)
117                  strdmp(ji,jj,jk) = zsa
118                  ttrdmp(ji,jj,jk) = zta
119               END DO
120            END DO
121         END DO
122         !
123      CASE ( 1 )                  !==  no damping in the turbocline (avt > 5 cm2/s)  ==!
124         DO jk = 1, jpkm1
125            DO jj = 2, jpjm1
126               DO ji = fs_2, fs_jpim1   ! vector opt.
127                  IF( avt(ji,jj,jk) <= 5.e-4 ) THEN
128                     zta = resto(ji,jj,jk) * ( t_dta(ji,jj,jk) - tsb(ji,jj,jk,jp_tem) )
129                     zsa = resto(ji,jj,jk) * ( s_dta(ji,jj,jk) - tsb(ji,jj,jk,jp_sal) )
130                  ELSE
131                     zta = 0.e0
132                     zsa = 0.e0   
133                  ENDIF
134                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta
135                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
136                  ! save the salinity trend (used in asmtrj)
137                  strdmp(ji,jj,jk) = zsa
138                  ttrdmp(ji,jj,jk) = zta
139               END DO
140            END DO
141         END DO
142         !
143      CASE ( 2 )                  !==  no damping in the mixed layer   ==!
144         DO jk = 1, jpkm1
145            DO jj = 2, jpjm1
146               DO ji = fs_2, fs_jpim1   ! vector opt.
147                  IF( fsdept(ji,jj,jk) >= hmlp (ji,jj) ) THEN
148                     zta = resto(ji,jj,jk) * ( t_dta(ji,jj,jk) - tsb(ji,jj,jk,jp_tem) )
149                     zsa = resto(ji,jj,jk) * ( s_dta(ji,jj,jk) - tsb(ji,jj,jk,jp_sal) )
150                  ELSE
151                     zta = 0.e0
152                     zsa = 0.e0   
153                  ENDIF
154                  ! add the trends to the general tracer trends
155                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + zta
156                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + zsa
157                  ! save the salinity trend (used in asmtrj)
158                  strdmp(ji,jj,jk) = zsa
159                  ttrdmp(ji,jj,jk) = zta
160               END DO
161            END DO
162         END DO
163         !
164      END SELECT
165
166      IF( l_trdtra )   THEN       ! trend diagnostic
167         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
168         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
169         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_dmp, ztrdt )
170         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_dmp, ztrds )
171         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds )
172      ENDIF
173      !                           ! Control print
174      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' dmp  - Ta: ', mask1=tmask,   &
175         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
176      !
177   END SUBROUTINE tra_dmp
178
179
180   SUBROUTINE tra_dmp_init
181      !!----------------------------------------------------------------------
182      !!                  ***  ROUTINE tra_dmp_init  ***
183      !!
184      !! ** Purpose :   Initialization for the newtonian damping
185      !!
186      !! ** Method  :   read the nammbf namelist and check the parameters
187      !!----------------------------------------------------------------------
188      NAMELIST/namtra_dmp/ nn_hdmp, nn_zdmp, rn_surf, rn_bot, rn_dep, nn_file
189      !!----------------------------------------------------------------------
190
191      REWIND ( numnam )                  ! Read Namelist namtra_dmp : temperature and salinity damping term
192      READ   ( numnam, namtra_dmp )
193      IF( lzoom )   nn_zdmp = 0           ! restoring to climatology at closed north or south boundaries
194
195      IF(lwp) THEN                       ! Namelist print
196         WRITE(numout,*)
197         WRITE(numout,*) 'tra_dmp : T and S newtonian damping'
198         WRITE(numout,*) '~~~~~~~'
199         WRITE(numout,*) '   Namelist namtra_dmp : set damping parameter'
200         WRITE(numout,*) '      T and S damping option         nn_hdmp = ', nn_hdmp
201         WRITE(numout,*) '      mixed layer damping option     nn_zdmp = ', nn_zdmp, '(zoom: forced to 0)'
202         WRITE(numout,*) '      surface time scale (days)      rn_surf = ', rn_surf
203         WRITE(numout,*) '      bottom time scale (days)       rn_bot  = ', rn_bot
204         WRITE(numout,*) '      depth of transition (meters)   rn_dep  = ', rn_dep
205         WRITE(numout,*) '      create a damping.coeff file    nn_file = ', nn_file
206      ENDIF
207
208      SELECT CASE ( nn_hdmp )
209      CASE (  -1  )   ;   IF(lwp) WRITE(numout,*) '   tracer damping in the Med & Red seas only'
210      CASE ( 1:90 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping poleward of', nn_hdmp, ' degrees'
211      CASE DEFAULT
212         WRITE(ctmp1,*) '          bad flag value for nn_hdmp = ', nn_hdmp
213         CALL ctl_stop(ctmp1)
214      END SELECT
215
216      SELECT CASE ( nn_zdmp )
217      CASE ( 0 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping throughout the water column'
218      CASE ( 1 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the turbocline (avt > 5 cm2/s)'
219      CASE ( 2 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the mixed layer'
220      CASE DEFAULT
221         WRITE(ctmp1,*) 'bad flag value for nn_zdmp = ', nn_zdmp
222         CALL ctl_stop(ctmp1)
223      END SELECT
224
225      IF( .NOT.lk_dtasal .OR. .NOT.lk_dtatem )   &
226         &   CALL ctl_stop( 'no temperature and/or salinity data define key_dtatem and key_dtasal' )
227
228      strdmp(:,:,:) = 0.e0       ! internal damping salinity trend (used in asmtrj)
229      ttrdmp(:,:,:) = 0.e0
230      !                          ! Damping coefficients initialization
231      IF( lzoom ) THEN   ;   CALL dtacof_zoom( resto )
232      ELSE               ;   CALL dtacof( nn_hdmp, rn_surf, rn_bot, rn_dep,  &
233                             &            nn_file, 'TRA'  , resto            )
234      ENDIF
235      !
236   END SUBROUTINE tra_dmp_init
237
238
239   SUBROUTINE dtacof_zoom( presto )
240      !!----------------------------------------------------------------------
241      !!                  ***  ROUTINE dtacof_zoom  ***
242      !!
243      !! ** Purpose :   Compute the damping coefficient for zoom domain
244      !!
245      !! ** Method  : - set along closed boundary due to zoom a damping over
246      !!                6 points with a max time scale of 5 days.
247      !!              - ORCA arctic/antarctic zoom: set the damping along
248      !!                south/north boundary over a latitude strip.
249      !!
250      !! ** Action  : - resto, the damping coeff. for T and S
251      !!----------------------------------------------------------------------
252      !! * Arguments
253      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)  ::  presto     !: restoring coeff. (s-1)
254      !
255      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
256      REAL(wp) ::   zlat, zlat0, zlat1, zlat2   ! temporary scalar
257      REAL(wp), DIMENSION(6)  ::   zfact        ! temporary workspace
258      !!----------------------------------------------------------------------
259
260      zfact(1) =  1.
261      zfact(2) =  1. 
262      zfact(3) = 11./12.
263      zfact(4) =  8./12.
264      zfact(5) =  4./12.
265      zfact(6) =  1./12.
266      zfact(:) = zfact(:) / ( 5. * rday )    ! 5 days max restoring time scale
267
268      presto(:,:,:) = 0.e0
269
270      ! damping along the forced closed boundary over 6 grid-points
271      DO jn = 1, 6
272         IF( lzoom_w )   presto( mi0(jn+jpizoom):mi1(jn+jpizoom), : , : )                    = zfact(jn)   ! west  closed
273         IF( lzoom_s )   presto( : , mj0(jn+jpjzoom):mj1(jn+jpjzoom), : )                    = zfact(jn)   ! south closed
274         IF( lzoom_e )   presto( mi0(jpiglo+jpizoom-1-jn):mi1(jpiglo+jpizoom-1-jn) , : , : ) = zfact(jn)   ! east  closed
275         IF( lzoom_n )   presto( : , mj0(jpjglo+jpjzoom-1-jn):mj1(jpjglo+jpjzoom-1-jn) , : ) = zfact(jn)   ! north closed
276      END DO
277
278      !                                           ! ====================================================
279      IF( lzoom_arct .AND. lzoom_anta ) THEN      !  ORCA configuration : arctic zoom or antarctic zoom
280         !                                        ! ====================================================
281         IF(lwp) WRITE(numout,*)
282         IF(lwp .AND. lzoom_arct ) WRITE(numout,*) '              dtacof_zoom : ORCA    Arctic zoom'
283         IF(lwp .AND. lzoom_arct ) WRITE(numout,*) '              dtacof_zoom : ORCA Antarctic zoom'
284         IF(lwp) WRITE(numout,*)
285         !
286         !                          ! Initialization :
287         presto(:,:,:) = 0.e0
288         zlat0 = 10.                     ! zlat0 : latitude strip where resto decreases
289         zlat1 = 30.                     ! zlat1 : resto = 1 before zlat1
290         zlat2 = zlat1 + zlat0           ! zlat2 : resto decreases from 1 to 0 between zlat1 and zlat2
291
292         DO jk = 2, jpkm1           ! Compute arrays resto ; value for internal damping : 5 days
293            DO jj = 1, jpj
294               DO ji = 1, jpi
295                  zlat = ABS( gphit(ji,jj) )
296                  IF( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN
297                     presto(ji,jj,jk) = 0.5 * ( 1./(5.*rday) ) * ( 1. - cos(rpi*(zlat2-zlat)/zlat0) ) 
298                  ELSEIF( zlat < zlat1 ) THEN
299                     presto(ji,jj,jk) = 1./(5.*rday)
300                  ENDIF
301               END DO
302            END DO
303         END DO
304         !
305      ENDIF
306      !                             ! Mask resto array
307      presto(:,:,:) = presto(:,:,:) * tmask(:,:,:)
308      !
309   END SUBROUTINE dtacof_zoom
310
311
312   SUBROUTINE dtacof( kn_hdmp, pn_surf, pn_bot, pn_dep,  &
313      &               kn_file, cdtype , presto           )
314      !!----------------------------------------------------------------------
315      !!                  ***  ROUTINE dtacof  ***
316      !!
317      !! ** Purpose :   Compute the damping coefficient
318      !!
319      !! ** Method  :   Arrays defining the damping are computed for each grid
320      !!                point for temperature and salinity (resto)
321      !!                Damping depends on distance to coast, depth and latitude
322      !!
323      !! ** Action  : - resto, the damping coeff. for T and S
324      !!----------------------------------------------------------------------
325      USE iom
326      USE ioipsl
327      !! * Arguments
328      INTEGER                         , INTENT(in   )  ::  kn_hdmp    !: damping option
329      REAL(wp)                        , INTENT(in   )  ::  pn_surf    !: surface time scale (days)
330      REAL(wp)                        , INTENT(in   )  ::  pn_bot     !: bottom time scale (days)
331      REAL(wp)                        , INTENT(in   )  ::  pn_dep     !: depth of transition (meters)
332      INTEGER                         , INTENT(in   )  ::  kn_file    !: save the damping coef on a file or not
333      CHARACTER(len=3)                , INTENT(in   )  ::  cdtype     ! =TRA or TRC (tracer indicator)
334      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)  ::  presto     !: restoring coeff. (s-1)
335      !
336      INTEGER ::   ji, jj, jk                   ! dummy loop indices
337      INTEGER ::   ii0, ii1, ij0, ij1           !    -          -
338      INTEGER ::   inum0                        ! logical unit for file restoring damping term
339      INTEGER ::   icot                         ! logical unit for file distance to the coast
340      REAL(wp) ::   zinfl, zlon                 ! temporary scalars
341      REAL(wp) ::   zlat, zlat0, zlat1, zlat2   !    -         -
342      REAL(wp) ::   zsdmp, zbdmp                !    -         -
343      REAL(wp), DIMENSION(jpk)         ::   zhfac
344      REAL(wp), DIMENSION(jpi,jpj)     ::   zmrs
345      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdct
346      CHARACTER(len=20)                ::   cfile
347      !!----------------------------------------------------------------------
348
349      ! ====================================
350      !  ORCA configuration : global domain
351      ! ====================================
352
353      IF(lwp) WRITE(numout,*)
354      IF(lwp) WRITE(numout,*) '              dtacof : Global domain of ORCA'
355      IF(lwp) WRITE(numout,*) '              ------------------------------'
356
357      ! ... Initialization :
358      presto(:,:,:) = 0.e0
359      !
360      IF( kn_hdmp > 0 ) THEN      !  Damping poleward of 'nn_hdmp' degrees  !
361         !                        !-----------------------------------------!
362         IF(lwp) WRITE(numout,*)
363         IF(lwp) WRITE(numout,*) '              Damping poleward of ', kn_hdmp,' deg.'
364         !
365         CALL iom_open ( 'dist.coast.nc', icot, ldstop = .FALSE. )
366         !
367         IF( icot > 0 ) THEN          ! distance-to-coast read in file
368            CALL iom_get  ( icot, jpdom_data, 'Tcoast', zdct )
369            CALL iom_close( icot )
370         ELSE                         ! distance-to-coast computed and saved in file (output in zdct)
371            CALL cofdis( zdct )
372         ENDIF
373
374         !                            ! Compute arrays resto
375         zinfl = 1000.e3                   ! distance of influence for damping term
376         zlat0 = 10.                       ! latitude strip where resto decreases
377         zlat1 = REAL( kn_hdmp )           ! resto = 0 between -zlat1 and zlat1
378         zlat2 = zlat1 + zlat0             ! resto increases from 0 to 1 between |zlat1| and |zlat2|
379
380         DO jj = 1, jpj
381            DO ji = 1, jpi
382               zlat = ABS( gphit(ji,jj) )
383               IF ( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN
384                  presto(ji,jj,1) = 0.5 * ( 1. - cos(rpi*(zlat-zlat1)/zlat0 ) )
385               ELSEIF ( zlat > zlat2 ) THEN
386                  presto(ji,jj,1) = 1.
387               ENDIF
388            END DO
389         END DO
390
391         IF ( kn_hdmp == 20 ) THEN       ! North Indian ocean (20N/30N x 45E/100E) : resto=0
392            DO jj = 1, jpj
393               DO ji = 1, jpi
394                  zlat = gphit(ji,jj)
395                  zlon = MOD( glamt(ji,jj), 360. )
396                  IF ( zlat1 < zlat .AND. zlat < zlat2 .AND. 45. < zlon .AND. zlon < 100. ) THEN
397                     presto(ji,jj,1) = 0.e0
398                  ENDIF
399               END DO
400            END DO
401         ENDIF
402
403         zsdmp = 1./(pn_surf * rday)
404         zbdmp = 1./(pn_bot  * rday)
405         DO jk = 2, jpkm1
406            DO jj = 1, jpj
407               DO ji = 1, jpi
408                  zdct(ji,jj,jk) = MIN( zinfl, zdct(ji,jj,jk) )
409                  !   ... Decrease the value in the vicinity of the coast
410                  presto(ji,jj,jk) = presto(ji,jj,1) * 0.5 * ( 1. - COS( rpi*zdct(ji,jj,jk)/zinfl) )
411                  !   ... Vertical variation from zsdmp (sea surface) to zbdmp (bottom)
412                  presto(ji,jj,jk) = presto(ji,jj,jk)      * ( zbdmp + (zsdmp-zbdmp)*EXP(-fsdept(ji,jj,jk)/pn_dep) )
413               END DO
414            END DO
415         END DO
416         !
417      ENDIF
418
419
420      IF( cp_cfg == "orca" .AND. ( kn_hdmp > 0 .OR. kn_hdmp == -1 ) ) THEN
421
422         !                                         ! =========================
423         !                                         !  Med and Red Sea damping
424         !                                         ! =========================
425         IF(lwp)WRITE(numout,*)
426         IF(lwp)WRITE(numout,*) '              ORCA configuration: Damping in Med and Red Seas'
427
428
429         zmrs(:,:) = 0.e0                             ! damping term on the Med or Red Sea
430
431         SELECT CASE ( jp_cfg )
432         !                                           ! =======================
433         CASE ( 4 )                                  !  ORCA_R4 configuration
434            !                                        ! =======================
435            ! Mediterranean Sea
436            ij0 =  50   ;   ij1 =  56
437            ii0 =  81   ;   ii1 =  91   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
438            ij0 =  50   ;   ij1 =  55
439            ii0 =  75   ;   ii1 =  80   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
440            ij0 =  52   ;   ij1 =  53
441            ii0 =  70   ;   ii1 =  74   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
442            ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea
443            DO jk = 1, 17
444               zhfac (jk) = 0.5*( 1.- COS( rpi*(jk-1)/16. ) ) / rday
445            END DO
446            DO jk = 18, jpkm1
447               zhfac (jk) = 1./rday
448            END DO
449            !                                        ! =======================
450         CASE ( 2 )                                  !  ORCA_R2 configuration
451            !                                        ! =======================
452            ! Mediterranean Sea
453            ij0 =  96   ;   ij1 = 110
454            ii0 = 157   ;   ii1 = 181   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
455            ij0 = 100   ;   ij1 = 110
456            ii0 = 144   ;   ii1 = 156   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
457            ij0 = 100   ;   ij1 = 103
458            ii0 = 139   ;   ii1 = 143   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
459            ! Decrease before Gibraltar Strait
460            ij0 = 101   ;   ij1 = 102
461            ii0 = 139   ;   ii1 = 141   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.e0
462            ii0 = 142   ;   ii1 = 142   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0 / 90.e0
463            ii0 = 143   ;   ii1 = 143   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40e0
464            ii0 = 144   ;   ii1 = 144   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75e0
465            ! Red Sea
466            ij0 =  87   ;   ij1 =  96
467            ii0 = 147   ;   ii1 = 163   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
468            ! Decrease before Bab el Mandeb Strait
469            ij0 =  91   ;   ij1 =  91
470            ii0 = 153   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.80e0
471            ij0 =  90   ;   ij1 =  90
472            ii0 = 153   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40e0
473            ij0 =  89   ;   ij1 =  89
474            ii0 = 158   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0 / 90.e0
475            ij0 =  88   ;   ij1 =  88
476            ii0 = 160   ;   ii1 = 163   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.e0
477            ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea
478            DO jk = 1, 17
479               zhfac (jk) = 0.5*( 1.- COS( rpi*(jk-1)/16. ) ) / rday
480            END DO
481            DO jk = 18, jpkm1
482               zhfac (jk) = 1./rday
483            END DO
484            !                                        ! =======================
485         CASE ( 05 )                                 !  ORCA_R05 configuration
486            !                                        ! =======================
487            ! Mediterranean Sea
488            ii0 = 568   ;   ii1 = 574 
489            ij0 = 324   ;   ij1 = 333   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
490            ii0 = 575   ;   ii1 = 658
491            ij0 = 314   ;   ij1 = 366   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
492            ! Black Sea (remaining part
493            ii0 = 641   ;   ii1 = 651
494            ij0 = 367   ;   ij1 = 372   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
495            ! Decrease before Gibraltar Strait
496            ij0 = 324   ;   ij1 = 333
497            ii0 = 565   ;   ii1 = 565   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0 / 90.e0
498            ii0 = 566   ;   ii1 = 566   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
499            ii0 = 567   ;   ii1 = 567   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75
500            ! Red Sea
501            ii0 = 641   ;   ii1 = 665
502            ij0 = 270   ;   ij1 = 310   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0
503            ! Decrease before Bab el Mandeb Strait
504            ii0 = 666   ;   ii1 = 675
505            ij0 = 270   ;   ij1 = 290   
506            DO ji = mi0(ii0), mi1(ii1)
507               zmrs( ji , mj0(ij0):mj1(ij1) ) = 0.1 * ABS( FLOAT(ji - mi1(ii1)) )
508            END DO
509            zsdmp = 1./(pn_surf * rday)
510            zbdmp = 1./(pn_bot  * rday)
511            DO jk = 1, jpk
512               zhfac (jk) = ( zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(1,1,jk)/pn_dep) )
513            END DO
514            !                                       ! ========================
515         CASE ( 025 )                               !  ORCA_R025 configuration
516            !                                       ! ========================
517            CALL ctl_stop( ' Not yet implemented in ORCA_R025' )
518            !
519         END SELECT
520
521         DO jk = 1, jpkm1
522            presto(:,:,jk) = zmrs(:,:) * zhfac(jk) + ( 1. - zmrs(:,:) ) * presto(:,:,jk)
523         END DO
524
525         ! Mask resto array and set to 0 first and last levels
526         presto(:,:, : ) = presto(:,:,:) * tmask(:,:,:)
527         presto(:,:, 1 ) = 0.e0
528         presto(:,:,jpk) = 0.e0
529         !                         !--------------------!
530      ELSE                         !     No damping     !
531         !                         !--------------------!
532         CALL ctl_stop( 'Choose a correct value of nn_hdmp or DO NOT defined key_tradmp' )
533      ENDIF
534
535      !                            !--------------------------------!
536      IF( kn_file == 1 ) THEN      !  save damping coef. in a file  !
537         !                         !--------------------------------!
538         IF(lwp) WRITE(numout,*) '              create damping.coeff.nc file'
539         IF( cdtype == 'TRA' ) cfile = 'damping.coeff'
540         IF( cdtype == 'TRC' ) cfile = 'damping.coeff.trc'
541         cfile = TRIM( cfile )
542         CALL iom_open  ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )
543         CALL iom_rstput( 0, 0, inum0, 'Resto', presto )
544         CALL iom_close ( inum0 )
545      ENDIF
546      !
547   END SUBROUTINE dtacof
548
549
550   SUBROUTINE cofdis( pdct )
551      !!----------------------------------------------------------------------
552      !!                 ***  ROUTINE cofdis  ***
553      !!
554      !! ** Purpose :   Compute the distance between ocean T-points and the
555      !!      ocean model coastlines. Save the distance in a NetCDF file.
556      !!
557      !! ** Method  :   For each model level, the distance-to-coast is
558      !!      computed as follows :
559      !!       - The coastline is defined as the serie of U-,V-,F-points
560      !!      that are at the ocean-land bound.
561      !!       - For each ocean T-point, the distance-to-coast is then
562      !!      computed as the smallest distance (on the sphere) between the
563      !!      T-point and all the coastline points.
564      !!       - For land T-points, the distance-to-coast is set to zero.
565      !!      C A U T I O N : Computation not yet implemented in mpp case.
566      !!
567      !! ** Action  : - pdct, distance to the coastline (argument)
568      !!              - NetCDF file 'dist.coast.nc'
569      !!----------------------------------------------------------------------
570      USE ioipsl      ! IOipsl librairy
571      !!
572      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) ::   pdct   ! distance to the coastline
573      !!
574      INTEGER ::   ji, jj, jk, jl      ! dummy loop indices
575      INTEGER ::   iju, ijt            ! temporary integers
576      INTEGER ::   icoast, itime
577      INTEGER ::   icot         ! logical unit for file distance to the coast
578      LOGICAL, DIMENSION(jpi,jpj) ::   llcotu, llcotv, llcotf   ! ???
579      CHARACTER (len=32) ::   clname
580      REAL(wp) ::   zdate0
581      REAL(wp), DIMENSION(jpi,jpj)   ::   zxt, zyt, zzt, zmask   ! cartesian coordinates for T-points
582      REAL(wp), DIMENSION(3*jpi*jpj) ::   zxc, zyc, zzc, zdis    ! temporary workspace
583      !!----------------------------------------------------------------------
584
585      ! 0. Initialization
586      ! -----------------
587      IF(lwp) WRITE(numout,*)
588      IF(lwp) WRITE(numout,*) 'cofdis : compute the distance to coastline'
589      IF(lwp) WRITE(numout,*) '~~~~~~'
590      IF(lwp) WRITE(numout,*)
591      IF( lk_mpp ) &
592           & CALL ctl_stop('         Computation not yet implemented with key_mpp_...', &
593           &               '         Rerun the code on another computer or ', &
594           &               '         create the "dist.coast.nc" file using IDL' )
595
596      pdct(:,:,:) = 0.e0
597      zxt(:,:) = cos( rad * gphit(:,:) ) * cos( rad * glamt(:,:) )
598      zyt(:,:) = cos( rad * gphit(:,:) ) * sin( rad * glamt(:,:) )
599      zzt(:,:) = sin( rad * gphit(:,:) )
600
601
602      ! 1. Loop on vertical levels
603      ! --------------------------
604      !                                                ! ===============
605      DO jk = 1, jpkm1                                 ! Horizontal slab
606         !                                             ! ===============
607         ! Define the coastline points (U, V and F)
608         DO jj = 2, jpjm1
609            DO ji = 2, jpim1
610               zmask(ji,jj) =  ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
611                   &           + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) )
612               llcotu(ji,jj) = ( tmask(ji,jj,  jk) + tmask(ji+1,jj  ,jk) == 1. ) 
613               llcotv(ji,jj) = ( tmask(ji,jj  ,jk) + tmask(ji  ,jj+1,jk) == 1. ) 
614               llcotf(ji,jj) = ( zmask(ji,jj) > 0. ) .AND. ( zmask(ji,jj) < 4. )
615            END DO
616         END DO
617
618         ! Lateral boundaries conditions
619         llcotu(:, 1 ) = umask(:,  2  ,jk) == 1
620         llcotu(:,jpj) = umask(:,jpjm1,jk) == 1
621         llcotv(:, 1 ) = vmask(:,  2  ,jk) == 1
622         llcotv(:,jpj) = vmask(:,jpjm1,jk) == 1
623         llcotf(:, 1 ) = fmask(:,  2  ,jk) == 1
624         llcotf(:,jpj) = fmask(:,jpjm1,jk) == 1
625
626         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
627            llcotu( 1 ,:) = llcotu(jpim1,:)
628            llcotu(jpi,:) = llcotu(  2  ,:)
629            llcotv( 1 ,:) = llcotv(jpim1,:)
630            llcotv(jpi,:) = llcotv(  2  ,:)
631            llcotf( 1 ,:) = llcotf(jpim1,:)
632            llcotf(jpi,:) = llcotf(  2  ,:)
633         ELSE
634            llcotu( 1 ,:) = umask(  2  ,:,jk) == 1
635            llcotu(jpi,:) = umask(jpim1,:,jk) == 1
636            llcotv( 1 ,:) = vmask(  2  ,:,jk) == 1
637            llcotv(jpi,:) = vmask(jpim1,:,jk) == 1
638            llcotf( 1 ,:) = fmask(  2  ,:,jk) == 1
639            llcotf(jpi,:) = fmask(jpim1,:,jk) == 1
640         ENDIF
641         IF( nperio == 3 .OR. nperio == 4 ) THEN
642            DO ji = 1, jpim1
643               iju = jpi - ji + 1
644               llcotu(ji,jpj  ) = llcotu(iju,jpj-2)
645               llcotf(ji,jpjm1) = llcotf(iju,jpj-2)
646               llcotf(ji,jpj  ) = llcotf(iju,jpj-3)
647            END DO
648            DO ji = jpi/2, jpim1
649               iju = jpi - ji + 1
650               llcotu(ji,jpjm1) = llcotu(iju,jpjm1)
651            END DO
652            DO ji = 2, jpi
653               ijt = jpi - ji + 2
654               llcotv(ji,jpjm1) = llcotv(ijt,jpj-2)
655               llcotv(ji,jpj  ) = llcotv(ijt,jpj-3)
656            END DO
657         ENDIF
658         IF( nperio == 5 .OR. nperio == 6 ) THEN
659            DO ji = 1, jpim1
660               iju = jpi - ji
661               llcotu(ji,jpj  ) = llcotu(iju,jpjm1)
662               llcotf(ji,jpj  ) = llcotf(iju,jpj-2)
663            END DO
664            DO ji = jpi/2, jpim1
665               iju = jpi - ji
666               llcotf(ji,jpjm1) = llcotf(iju,jpjm1)
667            END DO
668            DO ji = 1, jpi
669               ijt = jpi - ji + 1
670               llcotv(ji,jpj  ) = llcotv(ijt,jpjm1)
671            END DO
672            DO ji = jpi/2+1, jpi
673               ijt = jpi - ji + 1
674               llcotv(ji,jpjm1) = llcotv(ijt,jpjm1)
675            END DO
676         ENDIF
677
678         ! Compute cartesian coordinates of coastline points
679         ! and the number of coastline points
680
681         icoast = 0
682         DO jj = 1, jpj
683            DO ji = 1, jpi
684               IF( llcotf(ji,jj) ) THEN
685                  icoast = icoast + 1
686                  zxc(icoast) = COS( rad*gphif(ji,jj) ) * COS( rad*glamf(ji,jj) )
687                  zyc(icoast) = COS( rad*gphif(ji,jj) ) * SIN( rad*glamf(ji,jj) )
688                  zzc(icoast) = SIN( rad*gphif(ji,jj) )
689               ENDIF
690               IF( llcotu(ji,jj) ) THEN
691                  icoast = icoast+1
692                  zxc(icoast) = COS( rad*gphiu(ji,jj) ) * COS( rad*glamu(ji,jj) )
693                  zyc(icoast) = COS( rad*gphiu(ji,jj) ) * SIN( rad*glamu(ji,jj) )
694                  zzc(icoast) = SIN( rad*gphiu(ji,jj) )
695               ENDIF
696               IF( llcotv(ji,jj) ) THEN
697                  icoast = icoast+1
698                  zxc(icoast) = COS( rad*gphiv(ji,jj) ) * COS( rad*glamv(ji,jj) )
699                  zyc(icoast) = COS( rad*gphiv(ji,jj) ) * SIN( rad*glamv(ji,jj) )
700                  zzc(icoast) = SIN( rad*gphiv(ji,jj) )
701               ENDIF
702            END DO
703         END DO
704
705         ! Distance for the T-points
706
707         DO jj = 1, jpj
708            DO ji = 1, jpi
709               IF( tmask(ji,jj,jk) == 0. ) THEN
710                  pdct(ji,jj,jk) = 0.
711               ELSE
712                  DO jl = 1, icoast
713                     zdis(jl) = ( zxt(ji,jj) - zxc(jl) )**2   &
714                        &     + ( zyt(ji,jj) - zyc(jl) )**2   &
715                        &     + ( zzt(ji,jj) - zzc(jl) )**2
716                  END DO
717                  pdct(ji,jj,jk) = ra * SQRT( MINVAL( zdis(1:icoast) ) )
718               ENDIF
719            END DO
720         END DO
721         !                                                ! ===============
722      END DO                                              !   End of slab
723      !                                                   ! ===============
724
725
726      ! 2. Create the  distance to the coast file in NetCDF format
727      ! ----------------------------------------------------------   
728      clname = 'dist.coast'
729      itime = 0
730      CALL ymds2ju( 0     , 1      , 1     , 0.e0 , zdate0 )
731      CALL restini( 'NONE', jpi    , jpj   , glamt, gphit ,   &
732         &          jpk   , gdept_0, clname, itime, zdate0,   &
733         &          rdt   , icot                         )
734      CALL restput( icot, 'Tcoast', jpi, jpj, jpk, 0, pdct )
735      CALL restclo( icot )
736
737   END SUBROUTINE cofdis
738
739#else
740   !!----------------------------------------------------------------------
741   !!   Default key                                     NO internal damping
742   !!----------------------------------------------------------------------
743   LOGICAL , PUBLIC, PARAMETER ::   lk_tradmp = .FALSE.    !: internal damping flag
744CONTAINS
745   SUBROUTINE tra_dmp( kt )        ! Empty routine
746      WRITE(*,*) 'tra_dmp: You should not have seen this print! error?', kt
747   END SUBROUTINE tra_dmp
748   SUBROUTINE tra_dmp_init        ! Empty routine
749   END SUBROUTINE tra_dmp_init
750#endif
751
752   !!======================================================================
753END MODULE tradmp
Note: See TracBrowser for help on using the repository browser.