source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/TOP/TRP/trcdmp.F90 @ 11671

Last change on this file since 11671 was 11671, checked in by acc, 19 months ago

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.

  • Property svn:keywords set to Id
File size: 18.2 KB
Line 
1MODULE trcdmp
2   !!======================================================================
3   !!                       ***  MODULE  trcdmp  ***
4   !! Ocean physics: internal restoring trend on passive tracers
5   !!======================================================================
6   !! History :  OPA  !  1991-03  (O. Marti, G. Madec)  Original code
7   !!                 !  1996-01  (G. Madec) statement function for e3
8   !!                 !  1997-05  (H. Loukos)  adapted for passive tracers
9   !!    NEMO    9.0  !  2004-03  (C. Ethe)    free form + modules
10   !!            3.2  !  2007-02  (C. Deltel)  Diagnose ML trends for passive tracers
11   !!            3.3  !  2010-06  (C. Ethe, G. Madec) merge TRA-TRC
12   !!----------------------------------------------------------------------
13#if  defined key_top 
14   !!----------------------------------------------------------------------
15   !!   trc_dmp      : update the tracer trend with the internal damping
16   !!   trc_dmp_init : initialization, namlist read, parameters control
17   !!----------------------------------------------------------------------
18   USE oce_trc         ! ocean dynamics and tracers variables
19   USE trc             ! ocean passive tracers variables
20   USE trcdta
21   USE tradmp
22   USE trdtra
23   USE trd_oce
24   !
25   USE iom
26   USE prtctl_trc      ! Print control for debbuging
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC trc_dmp     
32   PUBLIC trc_dmp_clo   
33   PUBLIC trc_dmp_alloc 
34   PUBLIC trc_dmp_ini   
35
36   INTEGER            , PUBLIC ::   nn_zdmp_tr    !: = 0/1/2 flag for damping in the mixed layer
37   CHARACTER(LEN=200) , PUBLIC ::   cn_resto_tr   !: File containing restoration coefficient
38
39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   restotr   ! restoring coeff. on tracers (s-1)
40
41   INTEGER, PARAMETER         ::   npncts = 8       ! number of closed sea
42   INTEGER, DIMENSION(npncts) ::   nctsi1, nctsj1   ! south-west closed sea limits (i,j)
43   INTEGER, DIMENSION(npncts) ::   nctsi2, nctsj2   ! north-east closed sea limits (i,j)
44
45   !! * Substitutions
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
49   !! $Id$
50   !! Software governed by the CeCILL license (see ./LICENSE)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   INTEGER FUNCTION trc_dmp_alloc()
55      !!----------------------------------------------------------------------
56      !!                   ***  ROUTINE trc_dmp_alloc  ***
57      !!----------------------------------------------------------------------
58      ALLOCATE( restotr(jpi,jpj,jpk) , STAT=trc_dmp_alloc )
59      !
60      IF( trc_dmp_alloc /= 0 )   CALL ctl_warn('trc_dmp_alloc: failed to allocate array')
61      !
62   END FUNCTION trc_dmp_alloc
63
64
65   SUBROUTINE trc_dmp( kt )
66      !!----------------------------------------------------------------------
67      !!                   ***  ROUTINE trc_dmp  ***
68      !!                 
69      !! ** Purpose :   Compute the passive tracer trend due to a newtonian damping
70      !!      of the tracer field towards given data field and add it to the
71      !!      general tracer trends.
72      !!
73      !! ** Method  :   Newtonian damping towards trdta computed
74      !!      and add to the general tracer trends:
75      !!                     trn = tra + restotr * (trdta - trb)
76      !!         The trend is computed either throughout the water column
77      !!      (nlmdmptr=0) or in area of weak vertical mixing (nlmdmptr=1) or
78      !!      below the well mixed layer (nlmdmptr=2)
79      !!
80      !! ** Action  : - update the tracer trends tra with the newtonian
81      !!                damping trends.
82      !!              - save the trends ('key_trdmxl_trc')
83      !!----------------------------------------------------------------------
84      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
85      !
86      INTEGER ::   ji, jj, jk, jn, jl   ! dummy loop indices
87      CHARACTER (len=22) ::   charout
88      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrtrd
89      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrcdta   ! 3D  workspace
90      !!----------------------------------------------------------------------
91      !
92      IF( ln_timing )   CALL timing_start('trc_dmp')
93      !
94      IF( l_trdtrc )   ALLOCATE( ztrtrd(jpi,jpj,jpk) )   ! temporary save of trends
95      !
96      IF( nb_trcdta > 0 ) THEN  ! Initialisation of tracer from a file that may also be used for damping
97         !
98         ALLOCATE( ztrcdta(jpi,jpj,jpk) )    ! Memory allocation
99         !                                                          ! ===========
100         DO jn = 1, jptra                                           ! tracer loop
101            !                                                       ! ===========
102            IF( l_trdtrc ) ztrtrd(:,:,:) = tra(:,:,:,jn)    ! save trends
103            !
104            IF( ln_trc_ini(jn) ) THEN      ! update passive tracers arrays with input data read from file
105               !
106               jl = n_trc_index(jn) 
107               CALL trc_dta( kt, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000
108               !
109               SELECT CASE ( nn_zdmp_tr )
110               !
111               CASE( 0 )                !==  newtonian damping throughout the water column  ==!
112                  DO jk = 1, jpkm1
113                     DO jj = 2, jpjm1
114                        DO ji = fs_2, fs_jpim1   ! vector opt.
115                           tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + restotr(ji,jj,jk) * ( ztrcdta(ji,jj,jk) - trb(ji,jj,jk,jn) )
116                        END DO
117                     END DO
118                  END DO
119                  !
120               CASE ( 1 )                !==  no damping in the turbocline (avt > 5 cm2/s)  ==!
121                  DO jk = 1, jpkm1
122                     DO jj = 2, jpjm1
123                        DO ji = fs_2, fs_jpim1   ! vector opt.
124                           IF( avt(ji,jj,jk) <= avt_c )  THEN
125                              tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + restotr(ji,jj,jk) * ( ztrcdta(ji,jj,jk) - trb(ji,jj,jk,jn) )
126                           ENDIF
127                        END DO
128                     END DO
129                  END DO
130                  !
131               CASE ( 2 )               !==  no damping in the mixed layer   ==!
132                  DO jk = 1, jpkm1
133                     DO jj = 2, jpjm1
134                        DO ji = fs_2, fs_jpim1   ! vector opt.
135                           IF( gdept_n(ji,jj,jk) >= hmlp (ji,jj) ) THEN
136                              tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + restotr(ji,jj,jk) * ( ztrcdta(ji,jj,jk) - trb(ji,jj,jk,jn) )
137                           END IF
138                        END DO
139                     END DO
140                  END DO
141                 
142               END SELECT
143               !
144            ENDIF
145            !
146            IF( l_trdtrc ) THEN
147               ztrtrd(:,:,:) = tra(:,:,:,jn) -  ztrtrd(:,:,:)
148               CALL trd_tra( kt, 'TRC', jn, jptra_dmp, ztrtrd )
149            END IF
150            !                                                       ! ===========
151         END DO                                                     ! tracer loop
152         !                                                          ! ===========
153         DEALLOCATE( ztrcdta )
154      ENDIF
155      !
156      IF( l_trdtrc )  DEALLOCATE( ztrtrd )
157      !                                          ! print mean trends (used for debugging)
158      IF( ln_ctl ) THEN
159         WRITE(charout, FMT="('dmp ')")
160         CALL prt_ctl_trc_info(charout)
161         CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
162      ENDIF
163      !
164      IF( ln_timing )   CALL timing_stop('trc_dmp')
165      !
166   END SUBROUTINE trc_dmp
167
168
169   SUBROUTINE trc_dmp_ini
170      !!----------------------------------------------------------------------
171      !!                  ***  ROUTINE trc_dmp_ini  ***
172      !!
173      !! ** Purpose :   Initialization for the newtonian damping
174      !!
175      !! ** Method  :   read the nammbf namelist and check the parameters
176      !!              called by trc_dmp at the first timestep (nittrc000)
177      !!----------------------------------------------------------------------
178      INTEGER ::   ios, imask  ! local integers
179      !!
180      NAMELIST/namtrc_dmp/ nn_zdmp_tr , cn_resto_tr
181      !!----------------------------------------------------------------------
182      !
183      READ  ( numnat_ref, namtrc_dmp, IOSTAT = ios, ERR = 909)
184909   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtrc_dmp in reference namelist' )
185      READ  ( numnat_cfg, namtrc_dmp, IOSTAT = ios, ERR = 910)
186910   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtrc_dmp in configuration namelist' )
187      IF(lwm) WRITE ( numont, namtrc_dmp )
188
189      IF(lwp) THEN                       ! Namelist print
190         WRITE(numout,*)
191         WRITE(numout,*) 'trc_dmp : Passive tracers newtonian damping'
192         WRITE(numout,*) '~~~~~~~'
193         WRITE(numout,*) '   Namelist namtrc_dmp : set damping parameter'
194         WRITE(numout,*) '      mixed layer damping option     nn_zdmp_tr  = ', nn_zdmp_tr, '(zoom: forced to 0)'
195         WRITE(numout,*) '      Restoration coeff file         cn_resto_tr = ', cn_resto_tr
196      ENDIF
197      !                          ! Allocate arrays
198      IF( trc_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trc_dmp_ini: unable to allocate arrays' )
199      !
200      SELECT CASE ( nn_zdmp_tr )
201      CASE ( 0 )   ;   IF(lwp) WRITE(numout,*) '      ===>>   tracer damping throughout the water column'
202      CASE ( 1 )   ;   IF(lwp) WRITE(numout,*) '      ===>>   no tracer damping in the turbocline (avt > 5 cm2/s)'
203      CASE ( 2 )   ;   IF(lwp) WRITE(numout,*) '      ===>>   no tracer damping in the mixed layer'
204      CASE DEFAULT
205         WRITE(ctmp1,*) 'bad flag value for nn_zdmp_tr = ', nn_zdmp_tr
206         CALL ctl_stop(ctmp1)
207      END SELECT
208
209      IF( .NOT.lk_c1d ) THEN
210         IF( .NOT.ln_tradmp )   &
211            &   CALL ctl_stop( 'passive tracer damping need ln_tradmp to compute damping coef.' )
212         !
213         !                          ! Read damping coefficients from file
214         !Read in mask from file
215         CALL iom_open ( cn_resto_tr, imask)
216         CALL iom_get  ( imask, jpdom_autoglo, 'resto', restotr)
217         CALL iom_close( imask )
218         !
219      ENDIF
220      !
221   END SUBROUTINE trc_dmp_ini
222
223
224   SUBROUTINE trc_dmp_clo( kt )
225      !!---------------------------------------------------------------------
226      !!                  ***  ROUTINE trc_dmp_clo  ***
227      !!
228      !! ** Purpose :   Closed sea domain initialization
229      !!
230      !! ** Method  :   if a closed sea is located only in a model grid point
231      !!                we restore to initial data
232      !!
233      !! ** Action  :   nctsi1(), nctsj1() : south-west closed sea limits (i,j)
234      !!                nctsi2(), nctsj2() : north-east Closed sea limits (i,j)
235      !!----------------------------------------------------------------------
236      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
237      !
238      INTEGER :: ji , jj, jk, jn, jl, jc                    ! dummy loop indicesa
239      INTEGER :: isrow                                      ! local index
240      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrcdta       ! 3D  workspace
241      !!----------------------------------------------------------------------
242
243      IF( kt == nit000 ) THEN
244         ! initial values
245         nctsi1(:) = 1  ;  nctsi2(:) = 1
246         nctsj1(:) = 1  ;  nctsj2(:) = 1
247
248         ! set the closed seas (in data domain indices)
249         ! -------------------
250
251         IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA") THEN
252            !
253            SELECT CASE ( nn_cfg )
254            !                                           ! =======================
255            CASE ( 1 )                                  ! eORCA_R1 configuration
256            !                                           ! =======================
257            isrow = 332 - jpjglo
258            !
259            nctsi1(1)   = 333  ; nctsj1(1)   = 243 - isrow   ! Caspian Sea
260            nctsi2(1)   = 342  ; nctsj2(1)   = 274 - isrow
261            !                                       
262            nctsi1(2)   = 198  ; nctsj1(2)   = 258 - isrow   ! Lake Superior
263            nctsi2(2)   = 204  ; nctsj2(2)   = 262 - isrow
264            !                                         
265            nctsi1(3)   = 201  ; nctsj1(3)   = 250 - isrow   ! Lake Michigan
266            nctsi2(3)   = 203  ; nctsj2(3)   = 256 - isrow
267            !                                       
268            nctsi1(4)   = 204  ; nctsj1(4)   = 252 - isrow   ! Lake Huron
269            nctsi2(4)   = 209  ; nctsj2(4)   = 256 - isrow
270            !                                       
271            nctsi1(5)   = 206  ; nctsj1(5)   = 249 - isrow   ! Lake Erie
272            nctsi2(5)   = 209  ; nctsj2(5)   = 251 - isrow
273            !                                       
274            nctsi1(6)   = 210  ; nctsj1(6)   = 252 - isrow   ! Lake Ontario
275            nctsi2(6)   = 212  ; nctsj2(6)   = 252 - isrow
276            !                                       
277            nctsi1(7)   = 321  ; nctsj1(7)   = 180 - isrow   ! Victoria Lake
278            nctsi2(7)   = 322  ; nctsj2(7)   = 189 - isrow
279            !                                       
280            nctsi1(8)   = 297  ; nctsj1(8)   = 270 - isrow   ! Baltic Sea
281            nctsi2(8)   = 308  ; nctsj2(8)   = 293 - isrow
282            !                                       
283            !                                           ! =======================
284            CASE ( 2 )                                  !  ORCA_R2 configuration
285               !                                        ! =======================
286               !                                     
287               nctsi1(1)   =  11  ;  nctsj1(1)   = 103       ! Caspian Sea
288               nctsi2(1)   =  17  ;  nctsj2(1)   = 112
289               !                                     
290               nctsi1(2)   =  97  ;  nctsj1(2)   = 107       ! Great North American Lakes
291               nctsi2(2)   = 103  ;  nctsj2(2)   = 111
292               !                                     
293               nctsi1(3)   = 174  ;  nctsj1(3)   = 107       ! Black Sea 1 : west part of the Black Sea
294               nctsi2(3)   = 181  ;  nctsj2(3)   = 112
295              !                                     
296               nctsi1(4)   =   2  ;  nctsj1(4)   = 107      ! Black Sea 2 : est part of the Black Sea
297               nctsi2(4)   =   6  ;  nctsj2(4)   = 112
298               !                                     
299               nctsi1(5)   =  145 ;  nctsj1(5)   = 116       ! Baltic Sea
300               nctsi2(5)   =  150 ;  nctsj2(5)   = 126
301               !                                        ! =======================
302            CASE ( 4 )                                  !  ORCA_R4 configuration
303               !                                        ! =======================
304               !                                   
305               nctsi1(1)   =  4  ;  nctsj1(1)   = 53         ! Caspian Sea
306               nctsi2(1)   =  4  ;  nctsj2(1)   = 56
307               !                                   
308               nctsi1(2)   = 49  ;  nctsj1(2)   = 55         ! Great North American Lakes
309               nctsi2(2)   = 51  ;  nctsj2(2)   = 56
310               !                                   
311               nctsi1(3)   = 88  ;  nctsj1(3)   = 55         ! Black Sea
312               nctsi2(3)   = 91  ;  nctsj2(3)   = 56
313               !                                   
314               nctsi1(4)   = 75  ;  nctsj1(4)   = 59         ! Baltic Sea
315               nctsi2(4)   = 76  ;  nctsj2(4)   = 61
316               !                                        ! =======================
317            CASE ( 025 )                                ! ORCA_R025 configuration
318               !                                        ! =======================
319               !                                   
320               nctsi1(1)   = 1330 ; nctsj1(1)   = 645        ! Caspian + Aral sea
321               nctsi2(1)   = 1400 ; nctsj2(1)   = 795
322               !                                   
323               nctsi1(2)   = 1284 ; nctsj1(2)   = 722        ! Azov Sea
324               nctsi2(2)   = 1304 ; nctsj2(2)   = 747
325               !
326            END SELECT
327            !
328         ENDIF
329         !
330         ! convert the position in local domain indices
331         ! --------------------------------------------
332         DO jc = 1, npncts
333            nctsi1(jc)   = mi0( nctsi1(jc) )
334            nctsj1(jc)   = mj0( nctsj1(jc) )
335            !
336            nctsi2(jc)   = mi1( nctsi2(jc) )
337            nctsj2(jc)   = mj1( nctsj2(jc) )
338         END DO
339         !
340      ENDIF
341
342      ! Restore close seas values to initial data
343      IF( ln_trcdta .AND. nb_trcdta > 0 )  THEN   ! Initialisation of tracer from a file that may also be used for damping
344         !
345         IF(lwp)  WRITE(numout,*)
346         IF(lwp)  WRITE(numout,*) ' trc_dmp_clo : Restoring of nutrients on close seas at time-step kt = ', kt
347         IF(lwp)  WRITE(numout,*)
348         !
349         ALLOCATE( ztrcdta(jpi,jpj,jpk) )   ! Memory allocation
350         !
351         DO jn = 1, jptra
352            IF( ln_trc_ini(jn) ) THEN      ! update passive tracers arrays with input data read from file
353                jl = n_trc_index(jn)
354                CALL trc_dta( kt, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000
355                DO jc = 1, npncts
356                   DO jk = 1, jpkm1
357                      DO jj = nctsj1(jc), nctsj2(jc)
358                         DO ji = nctsi1(jc), nctsi2(jc)
359                            trn(ji,jj,jk,jn) = ztrcdta(ji,jj,jk)
360                            trb(ji,jj,jk,jn) = trn(ji,jj,jk,jn)
361                         END DO
362                      END DO
363                   END DO
364                END DO
365             ENDIF
366          END DO
367          DEALLOCATE( ztrcdta )
368      ENDIF
369      !
370   END SUBROUTINE trc_dmp_clo
371 
372#else
373   !!----------------------------------------------------------------------
374   !!  Dummy module :                                     No passive tracer
375   !!----------------------------------------------------------------------
376CONTAINS
377   SUBROUTINE trc_dmp( kt )        ! Empty routine
378      INTEGER, INTENT(in) :: kt
379      WRITE(*,*) 'trc_dmp: You should not have seen this print! error?', kt
380   END SUBROUTINE trc_dmp
381#endif
382
383   !!======================================================================
384END MODULE trcdmp
Note: See TracBrowser for help on using the repository browser.