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

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

Update NEMOGCM from branch nemo_v3_3_beta

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