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

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90 @ 4292

Last change on this file since 4292 was 4292, checked in by cetlod, 10 years ago

dev_MERGE_2013 : 1st step of the merge, see ticket #1185
  • Property svn:keywords set to Id
File size: 39.3 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   !: internal damping flag
52   INTEGER , PUBLIC ::   nn_hdmp     ! = 0/-1/'latitude' for damping over T and S
53   INTEGER , PUBLIC ::   nn_zdmp     ! = 0/1/2 flag for damping in the mixed layer
54   REAL(wp), PUBLIC ::   rn_surf     ! surface time scale for internal damping        [days]
55   REAL(wp), PUBLIC ::   rn_bot      ! bottom time scale for internal damping         [days]
56   REAL(wp), PUBLIC ::   rn_dep      ! depth of transition between rn_surf and rn_bot [meters]
57   INTEGER , PUBLIC ::   nn_file     ! = 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      INTEGER  ::   ios                 ! Local integer output status for namelist read
198      !!----------------------------------------------------------------------
199
200      REWIND( numnam_ref )              ! Namelist namtra_dmp in reference namelist : Temperature and salinity damping term
201      READ  ( numnam_ref, namtra_dmp, IOSTAT = ios, ERR = 901)
202901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in reference namelist', lwp )
203
204      REWIND( numnam_cfg )              ! Namelist namtra_dmp in configuration namelist : Temperature and salinity damping term
205      READ  ( numnam_cfg, namtra_dmp, IOSTAT = ios, ERR = 902 )
206902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist', lwp )
207      WRITE ( numond, namtra_dmp )
208     
209      IF( lzoom .AND. .NOT. lk_c1d )   nn_zdmp = 0          ! restoring to climatology at closed north or south boundaries
210
211      IF(lwp) THEN                       ! Namelist print
212         WRITE(numout,*)
213         WRITE(numout,*) 'tra_dmp_init : T and S newtonian damping'
214         WRITE(numout,*) '~~~~~~~'
215         WRITE(numout,*) '   Namelist namtra_dmp : set damping parameter'
216         WRITE(numout,*) '      add a damping term or not       ln_tradmp = ', ln_tradmp
217         WRITE(numout,*) '      T and S damping option          nn_hdmp   = ', nn_hdmp
218         WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp, '(non-C1D zoom: forced to 0)'
219         WRITE(numout,*) '      surface time scale (days)       rn_surf   = ', rn_surf
220         WRITE(numout,*) '      bottom time scale (days)        rn_bot    = ', rn_bot
221         WRITE(numout,*) '      depth of transition (meters)    rn_dep    = ', rn_dep
222         WRITE(numout,*) '      create a damping.coeff file     nn_file   = ', nn_file
223         WRITE(numout,*)
224      ENDIF
225
226      IF( ln_tradmp ) THEN               ! initialization for T-S damping
227         !
228         IF( tra_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_dmp_init: unable to allocate arrays' )
229         !
230#if ! defined key_c1d
231         SELECT CASE ( nn_hdmp )
232         CASE (  -1  )   ;   IF(lwp) WRITE(numout,*) '   tracer damping in the Med & Red seas only'
233         CASE ( 1:90 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping poleward of', nn_hdmp, ' degrees'
234         CASE DEFAULT
235            WRITE(ctmp1,*) '          bad flag value for nn_hdmp = ', nn_hdmp
236            CALL ctl_stop(ctmp1)
237         END SELECT
238         !
239#endif
240         SELECT CASE ( nn_zdmp )
241         CASE ( 0 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping throughout the water column'
242         CASE ( 1 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the turbocline (avt > 5 cm2/s)'
243         CASE ( 2 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the mixed layer'
244         CASE DEFAULT
245            WRITE(ctmp1,*) 'bad flag value for nn_zdmp = ', nn_zdmp
246            CALL ctl_stop(ctmp1)
247         END SELECT
248         !
249         IF( .NOT.ln_tsd_tradmp ) THEN
250            CALL ctl_warn( 'tra_dmp_init: read T-S data not initialized, we force ln_tsd_tradmp=T' )
251            CALL dta_tsd_init( ld_tradmp=ln_tradmp )        ! forces the initialisation of T-S data
252         ENDIF
253         !
254         strdmp(:,:,:) = 0._wp       ! internal damping salinity trend (used in asmtrj)
255         ttrdmp(:,:,:) = 0._wp
256         !                          ! Damping coefficients initialization
257         IF( lzoom .AND. .NOT. lk_c1d ) THEN   ;   CALL dtacof_zoom( resto )
258         ELSE               ;   CALL dtacof( nn_hdmp, rn_surf, rn_bot, rn_dep, nn_file, 'TRA', resto )
259         ENDIF
260         !
261      ENDIF
262      !
263   END SUBROUTINE tra_dmp_init
264
265
266   SUBROUTINE dtacof_zoom( presto )
267      !!----------------------------------------------------------------------
268      !!                  ***  ROUTINE dtacof_zoom  ***
269      !!
270      !! ** Purpose :   Compute the damping coefficient for zoom domain
271      !!
272      !! ** Method  : - set along closed boundary due to zoom a damping over
273      !!                6 points with a max time scale of 5 days.
274      !!              - ORCA arctic/antarctic zoom: set the damping along
275      !!                south/north boundary over a latitude strip.
276      !!
277      !! ** Action  : - resto, the damping coeff. for T and S
278      !!----------------------------------------------------------------------
279      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)  ::   presto   ! restoring coeff. (s-1)
280      !
281      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
282      REAL(wp) ::   zlat, zlat0, zlat1, zlat2, z1_5d   ! local scalar
283      REAL(wp), DIMENSION(6)  ::   zfact               ! 1Dworkspace
284      !!----------------------------------------------------------------------
285      !
286      IF( nn_timing == 1 )  CALL timing_start( 'dtacof_zoom')
287      !
288
289      zfact(1) =  1._wp
290      zfact(2) =  1._wp
291      zfact(3) = 11._wp / 12._wp
292      zfact(4) =  8._wp / 12._wp
293      zfact(5) =  4._wp / 12._wp
294      zfact(6) =  1._wp / 12._wp
295      zfact(:) = zfact(:) / ( 5._wp * rday )    ! 5 days max restoring time scale
296
297      presto(:,:,:) = 0._wp
298
299      ! damping along the forced closed boundary over 6 grid-points
300      DO jn = 1, 6
301         IF( lzoom_w )   presto( mi0(jn+jpizoom):mi1(jn+jpizoom), : , : )                    = zfact(jn)   ! west  closed
302         IF( lzoom_s )   presto( : , mj0(jn+jpjzoom):mj1(jn+jpjzoom), : )                    = zfact(jn)   ! south closed
303         IF( lzoom_e )   presto( mi0(jpiglo+jpizoom-1-jn):mi1(jpiglo+jpizoom-1-jn) , : , : ) = zfact(jn)   ! east  closed
304         IF( lzoom_n )   presto( : , mj0(jpjglo+jpjzoom-1-jn):mj1(jpjglo+jpjzoom-1-jn) , : ) = zfact(jn)   ! north closed
305      END DO
306
307      !                                           ! ====================================================
308      IF( cp_cfz == "arctic" .OR. cp_cfz == "antarctic" ) THEN   !  ORCA configuration : arctic or antarctic zoom
309         !                                        ! ====================================================
310         IF(lwp) WRITE(numout,*)
311         IF(lwp .AND. cp_cfz == "arctic" ) WRITE(numout,*) '              dtacof_zoom : ORCA    Arctic zoom'
312         IF(lwp .AND. cp_cfz == "antarctic" ) WRITE(numout,*) '           dtacof_zoom : ORCA Antarctic zoom'
313         IF(lwp) WRITE(numout,*)
314         !
315         !                          ! Initialization :
316         presto(:,:,:) = 0._wp
317         zlat0 = 10._wp                     ! zlat0 : latitude strip where resto decreases
318         zlat1 = 30._wp                     ! zlat1 : resto = 1 before zlat1
319         zlat2 = zlat1 + zlat0              ! zlat2 : resto decreases from 1 to 0 between zlat1 and zlat2
320         z1_5d = 1._wp / ( 5._wp * rday )   ! z1_5d : 1 / 5days
321
322         DO jk = 2, jpkm1           ! Compute arrays resto ; value for internal damping : 5 days
323            DO jj = 1, jpj
324               DO ji = 1, jpi
325                  zlat = ABS( gphit(ji,jj) )
326                  IF( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN
327                     presto(ji,jj,jk) = 0.5_wp * z1_5d * (  1._wp - COS( rpi*(zlat2-zlat)/zlat0 )  ) 
328                  ELSEIF( zlat < zlat1 ) THEN
329                     presto(ji,jj,jk) = z1_5d
330                  ENDIF
331               END DO
332            END DO
333         END DO
334         !
335      ENDIF
336      !                             ! Mask resto array
337      presto(:,:,:) = presto(:,:,:) * tmask(:,:,:)
338      !
339      IF( nn_timing == 1 )  CALL timing_stop( 'dtacof_zoom')
340      !
341   END SUBROUTINE dtacof_zoom
342
343
344   SUBROUTINE dtacof( kn_hdmp, pn_surf, pn_bot, pn_dep,  &
345      &               kn_file, cdtype , presto           )
346      !!----------------------------------------------------------------------
347      !!                  ***  ROUTINE dtacof  ***
348      !!
349      !! ** Purpose :   Compute the damping coefficient
350      !!
351      !! ** Method  :   Arrays defining the damping are computed for each grid
352      !!                point for temperature and salinity (resto)
353      !!                Damping depends on distance to coast, depth and latitude
354      !!
355      !! ** Action  : - resto, the damping coeff. for T and S
356      !!----------------------------------------------------------------------
357      USE iom
358      USE ioipsl
359      !!
360      INTEGER                         , INTENT(in   )  ::  kn_hdmp    ! damping option
361      REAL(wp)                        , INTENT(in   )  ::  pn_surf    ! surface time scale (days)
362      REAL(wp)                        , INTENT(in   )  ::  pn_bot     ! bottom time scale (days)
363      REAL(wp)                        , INTENT(in   )  ::  pn_dep     ! depth of transition (meters)
364      INTEGER                         , INTENT(in   )  ::  kn_file    ! save the damping coef on a file or not
365      CHARACTER(len=3)                , INTENT(in   )  ::  cdtype     ! =TRA, TRC or DYN (tracer/dynamics indicator)
366      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)  ::  presto     ! restoring coeff. (s-1)
367      !
368      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
369      INTEGER  ::   ii0, ii1, ij0, ij1          ! local integers
370      INTEGER  ::   inum0, icot                 !   -       -
371      REAL(wp) ::   zinfl, zlon                 ! local scalars
372      REAL(wp) ::   zlat, zlat0, zlat1, zlat2   !   -      -
373      REAL(wp) ::   zsdmp, zbdmp                !   -      -
374      CHARACTER(len=20)                   :: cfile
375      REAL(wp), POINTER, DIMENSION(:    ) :: zhfac 
376      REAL(wp), POINTER, DIMENSION(:,:  ) :: zmrs 
377      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdct 
378      !!----------------------------------------------------------------------
379      !
380      IF( nn_timing == 1 )  CALL timing_start('dtacof')
381      !
382      CALL wrk_alloc( jpk, zhfac          )
383      CALL wrk_alloc( jpi, jpj, zmrs      )
384      CALL wrk_alloc( jpi, jpj, jpk, zdct )
385#if defined key_c1d
386      !                                   ! ====================
387      !                                   !  C1D configuration : local domain
388      !                                   ! ====================
389      !
390      IF(lwp) WRITE(numout,*)
391      IF(lwp) WRITE(numout,*) '              dtacof : C1D 3x3 local domain'
392      IF(lwp) WRITE(numout,*) '              -----------------------------'
393      !
394      presto(:,:,:) = 0._wp
395      !
396      zsdmp = 1._wp / ( pn_surf * rday )
397      zbdmp = 1._wp / ( pn_bot  * rday )
398      DO jk = 2, jpkm1
399         DO jj = 1, jpj
400            DO ji = 1, jpi
401               !   ONLY vertical variation from zsdmp (sea surface) to zbdmp (bottom)
402               presto(ji,jj,jk) = zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(ji,jj,jk)/pn_dep)
403            END DO
404         END DO
405      END DO
406      !
407      presto(:,:, : ) = presto(:,:,:) * tmask(:,:,:)
408#else
409      !                                   ! ====================
410      !                                   !  ORCA configuration : global domain
411      !                                   ! ====================
412      !
413      IF(lwp) WRITE(numout,*)
414      IF(lwp) WRITE(numout,*) '              dtacof : Global domain of ORCA'
415      IF(lwp) WRITE(numout,*) '              ------------------------------'
416      !
417      presto(:,:,:) = 0._wp
418      !
419      IF( kn_hdmp > 0 ) THEN      !  Damping poleward of 'nn_hdmp' degrees  !
420         !                        !-----------------------------------------!
421         IF(lwp) WRITE(numout,*)
422         IF(lwp) WRITE(numout,*) '              Damping poleward of ', kn_hdmp, ' deg.'
423         !
424         CALL iom_open ( 'dist.coast.nc', icot, ldstop = .FALSE. )
425         !
426         IF( icot > 0 ) THEN          ! distance-to-coast read in file
427            CALL iom_get  ( icot, jpdom_data, 'Tcoast', zdct )
428            CALL iom_close( icot )
429         ELSE                         ! distance-to-coast computed and saved in file (output in zdct)
430            CALL cofdis( zdct )
431         ENDIF
432
433         !                            ! Compute arrays resto
434         zinfl = 1000.e3_wp                ! distance of influence for damping term
435         zlat0 = 10._wp                    ! latitude strip where resto decreases
436         zlat1 = REAL( kn_hdmp )           ! resto = 0 between -zlat1 and zlat1
437         zlat2 = zlat1 + zlat0             ! resto increases from 0 to 1 between |zlat1| and |zlat2|
438
439         DO jj = 1, jpj
440            DO ji = 1, jpi
441               zlat = ABS( gphit(ji,jj) )
442               IF ( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN
443                  presto(ji,jj,1) = 0.5_wp * (  1._wp - COS( rpi*(zlat-zlat1)/zlat0 )  )
444               ELSEIF ( zlat > zlat2 ) THEN
445                  presto(ji,jj,1) = 1._wp
446               ENDIF
447            END DO
448         END DO
449
450         IF ( kn_hdmp == 20 ) THEN       ! North Indian ocean (20N/30N x 45E/100E) : resto=0
451            DO jj = 1, jpj
452               DO ji = 1, jpi
453                  zlat = gphit(ji,jj)
454                  zlon = MOD( glamt(ji,jj), 360._wp )
455                  IF ( zlat1 < zlat .AND. zlat < zlat2 .AND. 45._wp < zlon .AND. zlon < 100._wp ) THEN
456                     presto(ji,jj,1) = 0._wp
457                  ENDIF
458               END DO
459            END DO
460         ENDIF
461
462         zsdmp = 1._wp / ( pn_surf * rday )
463         zbdmp = 1._wp / ( pn_bot  * rday )
464         DO jk = 2, jpkm1
465            DO jj = 1, jpj
466               DO ji = 1, jpi
467                  zdct(ji,jj,jk) = MIN( zinfl, zdct(ji,jj,jk) )
468                  !   ... Decrease the value in the vicinity of the coast
469                  presto(ji,jj,jk) = presto(ji,jj,1 ) * 0.5_wp * (  1._wp - COS( rpi*zdct(ji,jj,jk)/zinfl)  )
470                  !   ... Vertical variation from zsdmp (sea surface) to zbdmp (bottom)
471                  presto(ji,jj,jk) = presto(ji,jj,jk) * (  zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(ji,jj,jk)/pn_dep)  )
472               END DO
473            END DO
474         END DO
475         !
476      ENDIF
477
478      !                                  ! =========================
479      !                                  !  Med and Red Sea damping    (ORCA configuration only)
480      !                                  ! =========================
481      IF( cp_cfg == "orca" .AND. ( kn_hdmp > 0 .OR. kn_hdmp == -1 ) ) THEN
482         IF(lwp)WRITE(numout,*)
483         IF(lwp)WRITE(numout,*) '              ORCA configuration: Damping in Med and Red Seas'
484         !
485         zmrs(:,:) = 0._wp
486         !
487         SELECT CASE ( jp_cfg )
488         !                                           ! =======================
489         CASE ( 4 )                                  !  ORCA_R4 configuration
490            !                                        ! =======================
491            ij0 =  50   ;   ij1 =  56                    ! Mediterranean Sea
492
493            ii0 =  81   ;   ii1 =  91   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
494            ij0 =  50   ;   ij1 =  55
495            ii0 =  75   ;   ii1 =  80   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
496            ij0 =  52   ;   ij1 =  53
497            ii0 =  70   ;   ii1 =  74   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
498            ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea
499            DO jk = 1, 17
500               zhfac (jk) = 0.5_wp * (  1._wp - COS( rpi * REAL(jk-1,wp) / 16._wp )  ) / rday
501            END DO
502            DO jk = 18, jpkm1
503               zhfac (jk) = 1._wp / rday
504            END DO
505            !                                        ! =======================
506         CASE ( 2 )                                  !  ORCA_R2 configuration
507            !                                        ! =======================
508            ij0 =  96   ;   ij1 = 110                    ! Mediterranean Sea
509            ii0 = 157   ;   ii1 = 181   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
510            ij0 = 100   ;   ij1 = 110
511            ii0 = 144   ;   ii1 = 156   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
512            ij0 = 100   ;   ij1 = 103
513            ii0 = 139   ;   ii1 = 143   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
514            !
515            ij0 = 101   ;   ij1 = 102                    ! Decrease before Gibraltar Strait
516            ii0 = 139   ;   ii1 = 141   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0._wp
517            ii0 = 142   ;   ii1 = 142   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp
518            ii0 = 143   ;   ii1 = 143   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp
519            ii0 = 144   ;   ii1 = 144   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75_wp
520            !
521            ij0 =  87   ;   ij1 =  96                    ! Red Sea
522            ii0 = 147   ;   ii1 = 163   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
523            !
524            ij0 =  91   ;   ij1 =  91                    ! Decrease before Bab el Mandeb Strait
525            ii0 = 153   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.80_wp
526            ij0 =  90   ;   ij1 =  90
527            ii0 = 153   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp
528            ij0 =  89   ;   ij1 =  89
529            ii0 = 158   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp
530            ij0 =  88   ;   ij1 =  88
531            ii0 = 160   ;   ii1 = 163   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0._wp
532            ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea
533            DO jk = 1, 17
534               zhfac (jk) = 0.5_wp * (  1._wp - COS( rpi * REAL(jk-1,wp) / 16._wp )  ) / rday
535            END DO
536            DO jk = 18, jpkm1
537               zhfac (jk) = 1._wp / rday
538            END DO
539            !                                        ! =======================
540         CASE ( 05 )                                 !  ORCA_R05 configuration
541            !                                        ! =======================
542            ii0 = 568   ;   ii1 = 574                    ! Mediterranean Sea
543            ij0 = 324   ;   ij1 = 333   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
544            ii0 = 575   ;   ii1 = 658
545            ij0 = 314   ;   ij1 = 366   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
546            !
547            ii0 = 641   ;   ii1 = 651                    ! Black Sea (remaining part
548            ij0 = 367   ;   ij1 = 372   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
549            !
550            ij0 = 324   ;   ij1 = 333                    ! Decrease before Gibraltar Strait
551            ii0 = 565   ;   ii1 = 565   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp
552            ii0 = 566   ;   ii1 = 566   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp
553            ii0 = 567   ;   ii1 = 567   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75_wp
554            !
555            ii0 = 641   ;   ii1 = 665                    ! Red Sea
556            ij0 = 270   ;   ij1 = 310   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp
557            !
558            ii0 = 666   ;   ii1 = 675                    ! Decrease before Bab el Mandeb Strait
559            ij0 = 270   ;   ij1 = 290   
560            DO ji = mi0(ii0), mi1(ii1)
561               zmrs( ji , mj0(ij0):mj1(ij1) ) = 0.1_wp * ABS( FLOAT(ji - mi1(ii1)) )
562            END DO
563            zsdmp = 1._wp / ( pn_surf * rday )
564            zbdmp = 1._wp / ( pn_bot  * rday )
565            DO jk = 1, jpk
566               zhfac(jk) = (  zbdmp + (zsdmp-zbdmp) * EXP( -fsdept(1,1,jk)/pn_dep )  )
567            END DO
568            !                                       ! ========================
569         CASE ( 025 )                               !  ORCA_R025 configuration
570            !                                       ! ========================
571            CALL ctl_stop( ' Not yet implemented in ORCA_R025' )
572            !
573         END SELECT
574
575         DO jk = 1, jpkm1
576            presto(:,:,jk) = zmrs(:,:) * zhfac(jk) + ( 1._wp - zmrs(:,:) ) * presto(:,:,jk)
577         END DO
578
579         ! Mask resto array and set to 0 first and last levels
580         presto(:,:, : ) = presto(:,:,:) * tmask(:,:,:)
581         presto(:,:, 1 ) = 0._wp
582         presto(:,:,jpk) = 0._wp
583         !                         !--------------------!
584      ELSE                         !     No damping     !
585         !                         !--------------------!
586         CALL ctl_stop( 'Choose a correct value of nn_hdmp or put ln_tradmp to FALSE' )
587      ENDIF
588#endif
589
590      !                            !--------------------------------!
591      IF( kn_file == 1 ) THEN      !  save damping coef. in a file  !
592         !                         !--------------------------------!
593         IF(lwp) WRITE(numout,*) '              create damping.coeff.nc file'
594         IF( cdtype == 'TRA' ) cfile = 'damping.coeff'
595         IF( cdtype == 'TRC' ) cfile = 'damping.coeff.trc'
596         IF( cdtype == 'DYN' ) cfile = 'damping.coeff.dyn'
597         cfile = TRIM( cfile )
598         CALL iom_open  ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )
599         CALL iom_rstput( 0, 0, inum0, 'Resto', presto )
600         CALL iom_close ( inum0 )
601      ENDIF
602      !
603      CALL wrk_dealloc( jpk, zhfac)
604      CALL wrk_dealloc( jpi, jpj, zmrs )
605      CALL wrk_dealloc( jpi, jpj, jpk, zdct )
606      !
607      IF( nn_timing == 1 )  CALL timing_stop('dtacof')
608      !
609   END SUBROUTINE dtacof
610
611
612   SUBROUTINE cofdis( pdct )
613      !!----------------------------------------------------------------------
614      !!                 ***  ROUTINE cofdis  ***
615      !!
616      !! ** Purpose :   Compute the distance between ocean T-points and the
617      !!      ocean model coastlines. Save the distance in a NetCDF file.
618      !!
619      !! ** Method  :   For each model level, the distance-to-coast is
620      !!      computed as follows :
621      !!       - The coastline is defined as the serie of U-,V-,F-points
622      !!      that are at the ocean-land bound.
623      !!       - For each ocean T-point, the distance-to-coast is then
624      !!      computed as the smallest distance (on the sphere) between the
625      !!      T-point and all the coastline points.
626      !!       - For land T-points, the distance-to-coast is set to zero.
627      !!      C A U T I O N : Computation not yet implemented in mpp case.
628      !!
629      !! ** Action  : - pdct, distance to the coastline (argument)
630      !!              - NetCDF file 'dist.coast.nc'
631      !!----------------------------------------------------------------------
632      USE ioipsl      ! IOipsl librairy
633      !!
634      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) ::   pdct   ! distance to the coastline
635      !!
636      INTEGER ::   ji, jj, jk, jl   ! dummy loop indices
637      INTEGER ::   iju, ijt, icoast, itime, ierr, icot   ! local integers
638      CHARACTER (len=32) ::   clname                     ! local name
639      REAL(wp) ::   zdate0                               ! local scalar
640      REAL(wp), POINTER, DIMENSION(:,:) ::  zxt, zyt, zzt, zmask
641      REAL(wp), POINTER, DIMENSION(:  ) ::  zxc, zyc, zzc, zdis    ! temporary workspace
642      LOGICAL , ALLOCATABLE, DIMENSION(:,:) ::  llcotu, llcotv, llcotf   ! 2D logical workspace
643      !!----------------------------------------------------------------------
644      !
645      IF( nn_timing == 1 )  CALL timing_start('cofdis')
646      !
647      CALL wrk_alloc( jpi, jpj , zxt, zyt, zzt, zmask    )
648      CALL wrk_alloc( 3*jpi*jpj, zxc, zyc, zzc, zdis     )
649      ALLOCATE( llcotu(jpi,jpj), llcotv(jpi,jpj), llcotf(jpi,jpj)  )
650      !
651      IF( lk_mpp    )   CALL mpp_sum( ierr )
652      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'cofdis: requested local arrays unavailable')
653
654      ! 0. Initialization
655      ! -----------------
656      IF(lwp) WRITE(numout,*)
657      IF(lwp) WRITE(numout,*) 'cofdis : compute the distance to coastline'
658      IF(lwp) WRITE(numout,*) '~~~~~~'
659      IF(lwp) WRITE(numout,*)
660      IF( lk_mpp ) &
661           & CALL ctl_stop('         Computation not yet implemented with key_mpp_...', &
662           &               '         Rerun the code on another computer or ', &
663           &               '         create the "dist.coast.nc" file using IDL' )
664
665      pdct(:,:,:) = 0._wp
666      zxt(:,:) = COS( rad * gphit(:,:) ) * COS( rad * glamt(:,:) )
667      zyt(:,:) = COS( rad * gphit(:,:) ) * SIN( rad * glamt(:,:) )
668      zzt(:,:) = SIN( rad * gphit(:,:) )
669
670
671      ! 1. Loop on vertical levels
672      ! --------------------------
673      !                                                ! ===============
674      DO jk = 1, jpkm1                                 ! Horizontal slab
675         !                                             ! ===============
676         ! Define the coastline points (U, V and F)
677         DO jj = 2, jpjm1
678            DO ji = 2, jpim1
679               zmask(ji,jj) =  ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
680                   &           + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) )
681               llcotu(ji,jj) = ( tmask(ji,jj,  jk) + tmask(ji+1,jj  ,jk) == 1._wp ) 
682               llcotv(ji,jj) = ( tmask(ji,jj  ,jk) + tmask(ji  ,jj+1,jk) == 1._wp ) 
683               llcotf(ji,jj) = ( zmask(ji,jj) > 0._wp ) .AND. ( zmask(ji,jj) < 4._wp )
684            END DO
685         END DO
686
687         ! Lateral boundaries conditions
688         llcotu(:, 1 ) = umask(:,  2  ,jk) == 1
689         llcotu(:,jpj) = umask(:,jpjm1,jk) == 1
690         llcotv(:, 1 ) = vmask(:,  2  ,jk) == 1
691         llcotv(:,jpj) = vmask(:,jpjm1,jk) == 1
692         llcotf(:, 1 ) = fmask(:,  2  ,jk) == 1
693         llcotf(:,jpj) = fmask(:,jpjm1,jk) == 1
694
695         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
696            llcotu( 1 ,:) = llcotu(jpim1,:)
697            llcotu(jpi,:) = llcotu(  2  ,:)
698            llcotv( 1 ,:) = llcotv(jpim1,:)
699            llcotv(jpi,:) = llcotv(  2  ,:)
700            llcotf( 1 ,:) = llcotf(jpim1,:)
701            llcotf(jpi,:) = llcotf(  2  ,:)
702         ELSE
703            llcotu( 1 ,:) = umask(  2  ,:,jk) == 1
704            llcotu(jpi,:) = umask(jpim1,:,jk) == 1
705            llcotv( 1 ,:) = vmask(  2  ,:,jk) == 1
706            llcotv(jpi,:) = vmask(jpim1,:,jk) == 1
707            llcotf( 1 ,:) = fmask(  2  ,:,jk) == 1
708            llcotf(jpi,:) = fmask(jpim1,:,jk) == 1
709         ENDIF
710         IF( nperio == 3 .OR. nperio == 4 ) THEN
711            DO ji = 1, jpim1
712               iju = jpi - ji + 1
713               llcotu(ji,jpj  ) = llcotu(iju,jpj-2)
714               llcotf(ji,jpjm1) = llcotf(iju,jpj-2)
715               llcotf(ji,jpj  ) = llcotf(iju,jpj-3)
716            END DO
717            DO ji = jpi/2, jpim1
718               iju = jpi - ji + 1
719               llcotu(ji,jpjm1) = llcotu(iju,jpjm1)
720            END DO
721            DO ji = 2, jpi
722               ijt = jpi - ji + 2
723               llcotv(ji,jpjm1) = llcotv(ijt,jpj-2)
724               llcotv(ji,jpj  ) = llcotv(ijt,jpj-3)
725            END DO
726         ENDIF
727         IF( nperio == 5 .OR. nperio == 6 ) THEN
728            DO ji = 1, jpim1
729               iju = jpi - ji
730               llcotu(ji,jpj  ) = llcotu(iju,jpjm1)
731               llcotf(ji,jpj  ) = llcotf(iju,jpj-2)
732            END DO
733            DO ji = jpi/2, jpim1
734               iju = jpi - ji
735               llcotf(ji,jpjm1) = llcotf(iju,jpjm1)
736            END DO
737            DO ji = 1, jpi
738               ijt = jpi - ji + 1
739               llcotv(ji,jpj  ) = llcotv(ijt,jpjm1)
740            END DO
741            DO ji = jpi/2+1, jpi
742               ijt = jpi - ji + 1
743               llcotv(ji,jpjm1) = llcotv(ijt,jpjm1)
744            END DO
745         ENDIF
746
747         ! Compute cartesian coordinates of coastline points
748         ! and the number of coastline points
749         icoast = 0
750         DO jj = 1, jpj
751            DO ji = 1, jpi
752               IF( llcotf(ji,jj) ) THEN
753                  icoast = icoast + 1
754                  zxc(icoast) = COS( rad*gphif(ji,jj) ) * COS( rad*glamf(ji,jj) )
755                  zyc(icoast) = COS( rad*gphif(ji,jj) ) * SIN( rad*glamf(ji,jj) )
756                  zzc(icoast) = SIN( rad*gphif(ji,jj) )
757               ENDIF
758               IF( llcotu(ji,jj) ) THEN
759                  icoast = icoast+1
760                  zxc(icoast) = COS( rad*gphiu(ji,jj) ) * COS( rad*glamu(ji,jj) )
761                  zyc(icoast) = COS( rad*gphiu(ji,jj) ) * SIN( rad*glamu(ji,jj) )
762                  zzc(icoast) = SIN( rad*gphiu(ji,jj) )
763               ENDIF
764               IF( llcotv(ji,jj) ) THEN
765                  icoast = icoast+1
766                  zxc(icoast) = COS( rad*gphiv(ji,jj) ) * COS( rad*glamv(ji,jj) )
767                  zyc(icoast) = COS( rad*gphiv(ji,jj) ) * SIN( rad*glamv(ji,jj) )
768                  zzc(icoast) = SIN( rad*gphiv(ji,jj) )
769               ENDIF
770            END DO
771         END DO
772
773         ! Distance for the T-points
774         DO jj = 1, jpj
775            DO ji = 1, jpi
776               IF( tmask(ji,jj,jk) == 0._wp ) THEN
777                  pdct(ji,jj,jk) = 0._wp
778               ELSE
779                  DO jl = 1, icoast
780                     zdis(jl) = ( zxt(ji,jj) - zxc(jl) )**2   &
781                        &     + ( zyt(ji,jj) - zyc(jl) )**2   &
782                        &     + ( zzt(ji,jj) - zzc(jl) )**2
783                  END DO
784                  pdct(ji,jj,jk) = ra * SQRT( MINVAL( zdis(1:icoast) ) )
785               ENDIF
786            END DO
787         END DO
788         !                                                ! ===============
789      END DO                                              !   End of slab
790      !                                                   ! ===============
791
792
793      ! 2. Create the  distance to the coast file in NetCDF format
794      ! ----------------------------------------------------------   
795      clname = 'dist.coast'
796      itime  = 0
797      CALL ymds2ju( 0     , 1       , 1     , 0._wp , zdate0 )
798      CALL restini( 'NONE', jpi     , jpj   , glamt, gphit ,   &
799         &          jpk   , gdept_1d, clname, itime, zdate0,   &
800         &          rdt   , icot                         )
801      CALL restput( icot, 'Tcoast', jpi, jpj, jpk, 0, pdct )
802      CALL restclo( icot )
803      !
804      CALL wrk_dealloc( jpi, jpj , zxt, zyt, zzt, zmask    )
805      CALL wrk_dealloc( 3*jpi*jpj, zxc, zyc, zzc, zdis     )
806      DEALLOCATE( llcotu, llcotv, llcotf  )
807      !
808      IF( nn_timing == 1 )  CALL timing_stop('cofdis')
809      !
810   END SUBROUTINE cofdis
811   !!======================================================================
812END MODULE tradmp
Note: See TracBrowser for help on using the repository browser.