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

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90 @ 4147

Last change on this file since 4147 was 4147, checked in by cetlod, 10 years ago

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