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.
trcrad.F90 in NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/TOP/TRP – NEMO

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

Last change on this file since 11671 was 11671, checked in by acc, 5 years 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: 14.3 KB
Line 
1MODULE trcrad
2   !!======================================================================
3   !!                       ***  MODULE  trcrad  ***
4   !! Ocean passive tracers:  correction of negative concentrations
5   !!======================================================================
6   !! History :   -   !  01-01  (O. Aumont & E. Kestenare)  Original code
7   !!            1.0  !  04-03  (C. Ethe)  free form F90
8   !!----------------------------------------------------------------------
9#if defined key_top
10   !!----------------------------------------------------------------------
11   !!   'key_top'                                                TOP models
12   !!----------------------------------------------------------------------
13   !!   trc_rad    : correction of negative concentrations
14   !!----------------------------------------------------------------------
15   USE par_trc             ! need jptra, number of passive tracers
16   USE oce_trc             ! ocean dynamics and tracers variables
17   USE trc                 ! ocean passive tracers variables
18   USE trd_oce
19   USE trdtra
20   USE prtctl_trc          ! Print control for debbuging
21   USE lib_fortran
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC trc_rad     
27   PUBLIC trc_rad_ini 
28
29   LOGICAL , PUBLIC ::   ln_trcrad           !: flag to artificially correct negative concentrations
30   REAL(wp), DIMENSION(:,:), ALLOCATABLE::   gainmass
31
32   !!----------------------------------------------------------------------
33   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
34   !! $Id$
35   !! Software governed by the CeCILL license (see ./LICENSE)
36   !!----------------------------------------------------------------------
37CONTAINS
38
39   SUBROUTINE trc_rad( kt )
40      !!----------------------------------------------------------------------
41      !!                  ***  ROUTINE trc_rad  ***
42      !!
43      !! ** Purpose :   "crappy" routine to correct artificial negative
44      !!              concentrations due to isopycnal scheme
45      !!
46      !! ** Method  : - PISCES or LOBSTER: Set negative concentrations to zero
47      !!                while computing the corresponding tracer content that
48      !!                is added to the tracers. Then, adjust the tracer
49      !!                concentration using a multiplicative factor so that
50      !!                the total tracer concentration is preserved.
51      !!              - CFC: simply set to zero the negative CFC concentration
52      !!                (the total CFC content is not strictly preserved)
53      !!----------------------------------------------------------------------
54      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
55      !
56      CHARACTER (len=22) :: charout
57      !!----------------------------------------------------------------------
58      !
59      IF( ln_timing )   CALL timing_start('trc_rad')
60      !
61      IF( ln_age     )   CALL trc_rad_sms( kt, trb, trn, jp_age , jp_age                )  !  AGE
62      IF( ll_cfc     )   CALL trc_rad_sms( kt, trb, trn, jp_cfc0, jp_cfc1               )  !  CFC model
63      IF( ln_c14     )   CALL trc_rad_sms( kt, trb, trn, jp_c14 , jp_c14                )  !  C14
64      IF( ln_pisces  )   CALL trc_rad_sms( kt, trb, trn, jp_pcs0, jp_pcs1, cpreserv='Y' )  !  PISCES model
65      IF( ln_my_trc  )   CALL trc_rad_sms( kt, trb, trn, jp_myt0, jp_myt1               )  !  MY_TRC model
66      !
67      IF(ln_ctl) THEN      ! print mean trends (used for debugging)
68         WRITE(charout, FMT="('rad')")
69         CALL prt_ctl_trc_info( charout )
70         CALL prt_ctl_trc( tab4d=trn, mask=tmask, clinfo=ctrcnm )
71      ENDIF
72      !
73      IF( ln_timing )   CALL timing_stop('trc_rad')
74      !
75   END SUBROUTINE trc_rad
76
77
78   SUBROUTINE trc_rad_ini
79      !!---------------------------------------------------------------------
80      !!                  ***  ROUTINE trc _rad_ini ***
81      !!
82      !! ** Purpose :   read  namelist options
83      !!----------------------------------------------------------------------
84      INTEGER ::   ios   ! Local integer output status for namelist read
85      !!
86      NAMELIST/namtrc_rad/ ln_trcrad
87      !!----------------------------------------------------------------------
88      !
89      READ  ( numnat_ref, namtrc_rad, IOSTAT = ios, ERR = 907)
90907   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtrc_rad in reference namelist' )
91      READ  ( numnat_cfg, namtrc_rad, IOSTAT = ios, ERR = 908 )
92908   IF( ios > 0 )   CALL ctl_nam ( ios , 'namtrc_rad in configuration namelist' )
93      IF(lwm) WRITE( numont, namtrc_rad )
94
95      IF(lwp) THEN                     !   ! Control print
96         WRITE(numout,*)
97         WRITE(numout,*) 'trc_rad : Correct artificial negative concentrations '
98         WRITE(numout,*) '~~~~~~~ '
99         WRITE(numout,*) '   Namelist namtrc_rad : treatment of negative concentrations'
100         WRITE(numout,*) '      correct artificially negative concen. or not   ln_trcrad = ', ln_trcrad
101         WRITE(numout,*)
102         IF( ln_trcrad ) THEN   ;   WRITE(numout,*) '      ===>>   ensure the global tracer conservation'
103         ELSE                   ;   WRITE(numout,*) '      ===>>   NO strict global tracer conservation'     
104         ENDIF
105      ENDIF
106      !
107      ALLOCATE( gainmass(jptra,2) )
108      gainmass(:,:) = 0.
109      !
110   END SUBROUTINE trc_rad_ini
111
112
113   SUBROUTINE trc_rad_sms( kt, ptrb, ptrn, jp_sms0, jp_sms1, cpreserv )
114      !!-----------------------------------------------------------------------------
115      !!                  ***  ROUTINE trc_rad_sms  ***
116      !!
117      !! ** Purpose :   "crappy" routine to correct artificial negative
118      !!              concentrations due to isopycnal scheme
119      !!
120      !! ** Method  : 2 cases :
121      !!                - Set negative concentrations to zero while computing
122      !!                  the corresponding tracer content that is added to the
123      !!                  tracers. Then, adjust the tracer concentration using
124      !!                  a multiplicative factor so that the total tracer
125      !!                  concentration is preserved.
126      !!                - simply set to zero the negative CFC concentration
127      !!                  (the total content of concentration is not strictly preserved)
128      !!--------------------------------------------------------------------------------
129      INTEGER                                , INTENT(in   ) ::   kt                 ! ocean time-step index
130      INTEGER                                , INTENT(in   ) ::   jp_sms0, jp_sms1   ! First & last index of the passive tracer model
131      REAL(wp), DIMENSION (jpi,jpj,jpk,jptra), INTENT(inout) ::   ptrb    , ptrn     ! before and now traceur concentration
132      CHARACTER( len = 1), OPTIONAL          , INTENT(in   ) ::   cpreserv           ! flag to preserve content or not
133      !
134      INTEGER ::   ji, ji2, jj, jj2, jk, jn     ! dummy loop indices
135      INTEGER ::   icnt
136      LOGICAL ::   lldebug = .FALSE.            ! local logical
137      REAL(wp)::   zcoef, zs2rdt, ztotmass
138      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrneg, ztrpos
139      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrtrd   ! workspace arrays
140      !!----------------------------------------------------------------------
141      !
142      IF( l_trdtrc )   ALLOCATE( ztrtrd(jpi,jpj,jpk) )
143      zs2rdt = 1. / ( 2. * rdt * REAL( nn_dttrc, wp ) )
144      !
145      IF( PRESENT( cpreserv )  ) THEN     !==  total tracer concentration is preserved  ==!
146         !
147         ALLOCATE( ztrneg(1:jpi,1:jpj,jp_sms0:jp_sms1), ztrpos(1:jpi,1:jpj,jp_sms0:jp_sms1) )
148
149         DO jn = jp_sms0, jp_sms1
150            ztrneg(:,:,jn) = SUM( MIN( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:), dim = 3 )   ! sum of the negative values
151            ztrpos(:,:,jn) = SUM( MAX( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:), dim = 3 )   ! sum of the positive values
152         END DO
153         CALL sum3x3( ztrneg )
154         CALL sum3x3( ztrpos )
155         
156         DO jn = jp_sms0, jp_sms1
157            !
158            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrb(:,:,:,jn)                            ! save input trb for trend computation           
159            !
160            DO jk = 1, jpkm1
161               DO jj = 1, jpj
162                  DO ji = 1, jpi
163                     IF( ztrneg(ji,jj,jn) /= 0. ) THEN                                 ! if negative values over the 3x3 box
164                        !
165                        ptrb(ji,jj,jk,jn) = ptrb(ji,jj,jk,jn) * tmask(ji,jj,jk)   ! really needed?
166                        IF( ptrb(ji,jj,jk,jn) < 0. ) ptrb(ji,jj,jk,jn) = 0.       ! supress negative values
167                        IF( ptrb(ji,jj,jk,jn) > 0. ) THEN                         ! use positive values to compensate mass gain
168                           zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn)       ! ztrpos > 0 as ptrb > 0
169                           ptrb(ji,jj,jk,jn) = ptrb(ji,jj,jk,jn) * zcoef
170                           IF( zcoef < 0. ) THEN                                  ! if the compensation exceed the positive value
171                              gainmass(jn,1) = gainmass(jn,1) - ptrb(ji,jj,jk,jn) * cvol(ji,jj,jk)   ! we are adding mass...
172                              ptrb(ji,jj,jk,jn) = 0.                              ! limit the compensation to keep positive value
173                           ENDIF
174                        ENDIF
175                        !
176                     ENDIF
177                  END DO
178               END DO
179            END DO
180            !
181            IF( l_trdtrc ) THEN
182               ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
183               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrd )       ! Asselin-like trend handling
184            ENDIF
185            !
186         END DO
187 
188         IF( kt == nitend ) THEN
189            CALL mpp_sum( 'trcrad', gainmass(:,1) )
190            DO jn = jp_sms0, jp_sms1
191               IF( gainmass(jn,1) > 0. ) THEN
192                  ztotmass = glob_sum( 'trcrad', ptrb(:,:,:,jn) * cvol(:,:,:) )
193                  IF(lwp) WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrb, traceur ', jn  &
194                     &        , ' total mass : ', ztotmass, ', mass gain : ',  gainmass(jn,1)
195               END IF
196            END DO
197         ENDIF
198
199         DO jn = jp_sms0, jp_sms1
200            ztrneg(:,:,jn) = SUM( MIN( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:), dim = 3 )   ! sum of the negative values
201            ztrpos(:,:,jn) = SUM( MAX( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:), dim = 3 )   ! sum of the positive values
202         END DO
203         CALL sum3x3( ztrneg )
204         CALL sum3x3( ztrpos )
205         
206         DO jn = jp_sms0, jp_sms1
207            !
208            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrn(:,:,:,jn)                            ! save input trb for trend computation
209            !
210            DO jk = 1, jpkm1
211               DO jj = 1, jpj
212                  DO ji = 1, jpi
213                     IF( ztrneg(ji,jj,jn) /= 0. ) THEN                                 ! if negative values over the 3x3 box
214                        !
215                        ptrn(ji,jj,jk,jn) = ptrn(ji,jj,jk,jn) * tmask(ji,jj,jk)   ! really needed?
216                        IF( ptrn(ji,jj,jk,jn) < 0. ) ptrn(ji,jj,jk,jn) = 0.       ! supress negative values
217                        IF( ptrn(ji,jj,jk,jn) > 0. ) THEN                         ! use positive values to compensate mass gain
218                           zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn)       ! ztrpos > 0 as ptrb > 0
219                           ptrn(ji,jj,jk,jn) = ptrn(ji,jj,jk,jn) * zcoef
220                           IF( zcoef < 0. ) THEN                                  ! if the compensation exceed the positive value
221                              gainmass(jn,2) = gainmass(jn,2) - ptrn(ji,jj,jk,jn) * cvol(ji,jj,jk)   ! we are adding mass...
222                              ptrn(ji,jj,jk,jn) = 0.                              ! limit the compensation to keep positive value
223                           ENDIF
224                        ENDIF
225                        !
226                     ENDIF
227                  END DO
228               END DO
229            END DO
230            !
231            IF( l_trdtrc ) THEN
232               ztrtrd(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
233               CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrd )       ! standard     trend handling
234            ENDIF
235            !
236         END DO
237 
238         IF( kt == nitend ) THEN
239            CALL mpp_sum( 'trcrad', gainmass(:,2) )
240            DO jn = jp_sms0, jp_sms1
241               IF( gainmass(jn,2) > 0. ) THEN
242                  ztotmass = glob_sum( 'trcrad', ptrn(:,:,:,jn) * cvol(:,:,:) )
243                  WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrn, traceur ', jn  &
244                     &        , ' total mass : ', ztotmass, ', mass gain : ',  gainmass(jn,1)
245               END IF
246            END DO
247         ENDIF
248
249         DEALLOCATE( ztrneg, ztrpos )
250         !
251      ELSE                                !==  total CFC content is NOT strictly preserved  ==!
252         !
253         DO jn = jp_sms0, jp_sms1 
254            !
255            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrb(:,:,:,jn)                        ! save input trb for trend computation
256            !
257            WHERE( ptrb(:,:,:,jn) < 0. )   ptrb(:,:,:,jn) = 0.
258            !
259            IF( l_trdtrc ) THEN
260               ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
261               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrd )       ! Asselin-like trend handling
262            ENDIF
263            !
264            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrn(:,:,:,jn)                        ! save input trn for trend computation
265            !
266            WHERE( ptrn(:,:,:,jn) < 0. )   ptrn(:,:,:,jn) = 0.
267            !
268            IF( l_trdtrc ) THEN
269               ztrtrd(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
270               CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrd )       ! standard     trend handling
271            ENDIF
272            !
273         END DO
274         !
275      ENDIF
276      !
277      IF( l_trdtrc )  DEALLOCATE( ztrtrd )
278      !
279   END SUBROUTINE trc_rad_sms
280
281#else
282   !!----------------------------------------------------------------------
283   !!   Dummy module :                                         NO TOP model
284   !!----------------------------------------------------------------------
285CONTAINS
286   SUBROUTINE trc_rad( kt )              ! Empty routine
287      INTEGER, INTENT(in) ::   kt
288      WRITE(*,*) 'trc_rad: You should not have seen this print! error?', kt
289   END SUBROUTINE trc_rad
290#endif
291   
292   !!======================================================================
293END MODULE trcrad
Note: See TracBrowser for help on using the repository browser.