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

source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90 @ 3787

Last change on this file since 3787 was 3294, checked in by rblod, 12 years ago

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