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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90 @ 2690

Last change on this file since 2690 was 2690, checked in by gm, 13 years ago

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