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

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90 @ 3317

Last change on this file since 3317 was 3317, checked in by gm, 12 years ago

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