source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90 @ 4409

Last change on this file since 4409 was 4409, checked in by trackstand2, 7 years ago

Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.

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