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

source: tags/nemo_v3_2/nemo_v3_2/NEMO/OPA_SRC/TRA/tradmp.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

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