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

source: branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90 @ 5075

Last change on this file since 5075 was 5075, checked in by timgraham, 9 years ago

Upgraded branch to current head of trunk (r5072) so it can be used with the trunk

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