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

source: branches/2013/dev_r3987_UKMO6_C1D/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90 @ 4144

Last change on this file since 4144 was 4144, checked in by rfurner, 10 years ago

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