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

source: branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90 @ 2789

Last change on this file since 2789 was 2789, checked in by cetlod, 13 years ago

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