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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcrad.F90 @ 10946

Last change on this file since 10946 was 10946, checked in by acc, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert STO, TRD and USR modules and all knock on effects of these conversions. Note change to USR module may have implications for the TEST CASES (not tested yet). Standard SETTE tested only
  • Property svn:keywords set to Id
File size: 14.7 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, Kmm, Krhs )
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      INTEGER, INTENT(in) ::   Kmm, Krhs  ! time level indices
56      !
57      CHARACTER (len=22) :: charout
58      !!----------------------------------------------------------------------
59      !
60      IF( ln_timing )   CALL timing_start('trc_rad')
61      !
62      IF( ln_age     )   CALL trc_rad_sms( kt, Kmm, Krhs, trb, trn, jp_age , jp_age                )  !  AGE
63      IF( ll_cfc     )   CALL trc_rad_sms( kt, Kmm, Krhs, trb, trn, jp_cfc0, jp_cfc1               )  !  CFC model
64      IF( ln_c14     )   CALL trc_rad_sms( kt, Kmm, Krhs, trb, trn, jp_c14 , jp_c14                )  !  C14
65      IF( ln_pisces  )   CALL trc_rad_sms( kt, Kmm, Krhs, trb, trn, jp_pcs0, jp_pcs1, cpreserv='Y' )  !  PISCES model
66      IF( ln_my_trc  )   CALL trc_rad_sms( kt, Kmm, Krhs, trb, trn, jp_myt0, jp_myt1               )  !  MY_TRC model
67      !
68      IF(ln_ctl) THEN      ! print mean trends (used for debugging)
69         WRITE(charout, FMT="('rad')")
70         CALL prt_ctl_trc_info( charout )
71         CALL prt_ctl_trc( tab4d=trn, mask=tmask, clinfo=ctrcnm )
72      ENDIF
73      !
74      IF( ln_timing )   CALL timing_stop('trc_rad')
75      !
76   END SUBROUTINE trc_rad
77
78
79   SUBROUTINE trc_rad_ini
80      !!---------------------------------------------------------------------
81      !!                  ***  ROUTINE trc _rad_ini ***
82      !!
83      !! ** Purpose :   read  namelist options
84      !!----------------------------------------------------------------------
85      INTEGER ::   ios   ! Local integer output status for namelist read
86      !!
87      NAMELIST/namtrc_rad/ ln_trcrad
88      !!----------------------------------------------------------------------
89      !
90      REWIND( numnat_ref )              ! namtrc_rad in reference namelist
91      READ  ( numnat_ref, namtrc_rad, IOSTAT = ios, ERR = 907)
92907   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtrc_rad in reference namelist', lwp )
93      REWIND( numnat_cfg )              ! namtrc_rad in configuration namelist
94      READ  ( numnat_cfg, namtrc_rad, IOSTAT = ios, ERR = 908 )
95908   IF( ios > 0 )   CALL ctl_nam ( ios , 'namtrc_rad in configuration namelist', lwp )
96      IF(lwm) WRITE( numont, namtrc_rad )
97
98      IF(lwp) THEN                     !   ! Control print
99         WRITE(numout,*)
100         WRITE(numout,*) 'trc_rad : Correct artificial negative concentrations '
101         WRITE(numout,*) '~~~~~~~ '
102         WRITE(numout,*) '   Namelist namtrc_rad : treatment of negative concentrations'
103         WRITE(numout,*) '      correct artificially negative concen. or not   ln_trcrad = ', ln_trcrad
104         WRITE(numout,*)
105         IF( ln_trcrad ) THEN   ;   WRITE(numout,*) '      ===>>   ensure the global tracer conservation'
106         ELSE                   ;   WRITE(numout,*) '      ===>>   NO strict global tracer conservation'     
107         ENDIF
108      ENDIF
109      !
110      ALLOCATE( gainmass(jptra,2) )
111      gainmass(:,:) = 0.
112      !
113   END SUBROUTINE trc_rad_ini
114
115
116   SUBROUTINE trc_rad_sms( kt, Kmm, Krhs, ptrb, ptrn, jp_sms0, jp_sms1, cpreserv )
117      !!-----------------------------------------------------------------------------
118      !!                  ***  ROUTINE trc_rad_sms  ***
119      !!
120      !! ** Purpose :   "crappy" routine to correct artificial negative
121      !!              concentrations due to isopycnal scheme
122      !!
123      !! ** Method  : 2 cases :
124      !!                - Set negative concentrations to zero while computing
125      !!                  the corresponding tracer content that is added to the
126      !!                  tracers. Then, adjust the tracer concentration using
127      !!                  a multiplicative factor so that the total tracer
128      !!                  concentration is preserved.
129      !!                - simply set to zero the negative CFC concentration
130      !!                  (the total content of concentration is not strictly preserved)
131      !!--------------------------------------------------------------------------------
132      INTEGER                                , INTENT(in   ) ::   kt                 ! ocean time-step index
133      INTEGER                                , INTENT(in   ) ::   Kmm, Krhs          ! time level indices
134      INTEGER                                , INTENT(in   ) ::   jp_sms0, jp_sms1   ! First & last index of the passive tracer model
135      REAL(wp), DIMENSION (jpi,jpj,jpk,jptra), INTENT(inout) ::   ptrb    , ptrn     ! before and now traceur concentration
136      CHARACTER( len = 1), OPTIONAL          , INTENT(in   ) ::   cpreserv           ! flag to preserve content or not
137      !
138      INTEGER ::   ji, ji2, jj, jj2, jk, jn     ! dummy loop indices
139      INTEGER ::   icnt
140      LOGICAL ::   lldebug = .FALSE.            ! local logical
141      REAL(wp)::   zcoef, zs2rdt, ztotmass
142      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrneg, ztrpos
143      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrtrd   ! workspace arrays
144      !!----------------------------------------------------------------------
145      !
146      IF( l_trdtrc )   ALLOCATE( ztrtrd(jpi,jpj,jpk) )
147      zs2rdt = 1. / ( 2. * rdt * REAL( nn_dttrc, wp ) )
148      !
149      IF( PRESENT( cpreserv )  ) THEN     !==  total tracer concentration is preserved  ==!
150         !
151         ALLOCATE( ztrneg(1:jpi,1:jpj,jp_sms0:jp_sms1), ztrpos(1:jpi,1:jpj,jp_sms0:jp_sms1) )
152
153         DO jn = jp_sms0, jp_sms1
154            ztrneg(:,:,jn) = SUM( MIN( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:), dim = 3 )   ! sum of the negative values
155            ztrpos(:,:,jn) = SUM( MAX( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:), dim = 3 )   ! sum of the positive values
156         END DO
157         CALL sum3x3( ztrneg )
158         CALL sum3x3( ztrpos )
159         
160         DO jn = jp_sms0, jp_sms1
161            !
162            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrb(:,:,:,jn)                            ! save input trb for trend computation           
163            !
164            DO jk = 1, jpkm1
165               DO jj = 1, jpj
166                  DO ji = 1, jpi
167                     IF( ztrneg(ji,jj,jn) /= 0. ) THEN                                 ! if negative values over the 3x3 box
168                        !
169                        ptrb(ji,jj,jk,jn) = ptrb(ji,jj,jk,jn) * tmask(ji,jj,jk)   ! really needed?
170                        IF( ptrb(ji,jj,jk,jn) < 0. ) ptrb(ji,jj,jk,jn) = 0.       ! supress negative values
171                        IF( ptrb(ji,jj,jk,jn) > 0. ) THEN                         ! use positive values to compensate mass gain
172                           zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn)       ! ztrpos > 0 as ptrb > 0
173                           ptrb(ji,jj,jk,jn) = ptrb(ji,jj,jk,jn) * zcoef
174                           IF( zcoef < 0. ) THEN                                  ! if the compensation exceed the positive value
175                              gainmass(jn,1) = gainmass(jn,1) - ptrb(ji,jj,jk,jn) * cvol(ji,jj,jk)   ! we are adding mass...
176                              ptrb(ji,jj,jk,jn) = 0.                              ! limit the compensation to keep positive value
177                           ENDIF
178                        ENDIF
179                        !
180                     ENDIF
181                  END DO
182               END DO
183            END DO
184            !
185            IF( l_trdtrc ) THEN
186               ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
187               CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radb, ztrtrd )       ! Asselin-like trend handling
188            ENDIF
189            !
190         END DO
191 
192         IF( kt == nitend ) THEN
193            CALL mpp_sum( 'trcrad', gainmass(:,1) )
194            DO jn = jp_sms0, jp_sms1
195               IF( gainmass(jn,1) > 0. ) THEN
196                  ztotmass = glob_sum( 'trcrad', ptrb(:,:,:,jn) * cvol(:,:,:) )
197                  IF(lwp) WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrb, traceur ', jn  &
198                     &        , ' total mass : ', ztotmass, ', mass gain : ',  gainmass(jn,1)
199               END IF
200            END DO
201         ENDIF
202
203         DO jn = jp_sms0, jp_sms1
204            ztrneg(:,:,jn) = SUM( MIN( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:), dim = 3 )   ! sum of the negative values
205            ztrpos(:,:,jn) = SUM( MAX( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:), dim = 3 )   ! sum of the positive values
206         END DO
207         CALL sum3x3( ztrneg )
208         CALL sum3x3( ztrpos )
209         
210         DO jn = jp_sms0, jp_sms1
211            !
212            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrn(:,:,:,jn)                            ! save input trb for trend computation
213            !
214            DO jk = 1, jpkm1
215               DO jj = 1, jpj
216                  DO ji = 1, jpi
217                     IF( ztrneg(ji,jj,jn) /= 0. ) THEN                                 ! if negative values over the 3x3 box
218                        !
219                        ptrn(ji,jj,jk,jn) = ptrn(ji,jj,jk,jn) * tmask(ji,jj,jk)   ! really needed?
220                        IF( ptrn(ji,jj,jk,jn) < 0. ) ptrn(ji,jj,jk,jn) = 0.       ! supress negative values
221                        IF( ptrn(ji,jj,jk,jn) > 0. ) THEN                         ! use positive values to compensate mass gain
222                           zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn)       ! ztrpos > 0 as ptrb > 0
223                           ptrn(ji,jj,jk,jn) = ptrn(ji,jj,jk,jn) * zcoef
224                           IF( zcoef < 0. ) THEN                                  ! if the compensation exceed the positive value
225                              gainmass(jn,2) = gainmass(jn,2) - ptrn(ji,jj,jk,jn) * cvol(ji,jj,jk)   ! we are adding mass...
226                              ptrn(ji,jj,jk,jn) = 0.                              ! limit the compensation to keep positive value
227                           ENDIF
228                        ENDIF
229                        !
230                     ENDIF
231                  END DO
232               END DO
233            END DO
234            !
235            IF( l_trdtrc ) THEN
236               ztrtrd(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
237               CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radn, ztrtrd )       ! standard     trend handling
238            ENDIF
239            !
240         END DO
241 
242         IF( kt == nitend ) THEN
243            CALL mpp_sum( 'trcrad', gainmass(:,2) )
244            DO jn = jp_sms0, jp_sms1
245               IF( gainmass(jn,2) > 0. ) THEN
246                  ztotmass = glob_sum( 'trcrad', ptrn(:,:,:,jn) * cvol(:,:,:) )
247                  WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrn, traceur ', jn  &
248                     &        , ' total mass : ', ztotmass, ', mass gain : ',  gainmass(jn,1)
249               END IF
250            END DO
251         ENDIF
252
253         DEALLOCATE( ztrneg, ztrpos )
254         !
255      ELSE                                !==  total CFC content is NOT strictly preserved  ==!
256         !
257         DO jn = jp_sms0, jp_sms1 
258            !
259            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrb(:,:,:,jn)                        ! save input trb for trend computation
260            !
261            WHERE( ptrb(:,:,:,jn) < 0. )   ptrb(:,:,:,jn) = 0.
262            !
263            IF( l_trdtrc ) THEN
264               ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
265               CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radb, ztrtrd )       ! Asselin-like trend handling
266            ENDIF
267            !
268            IF( l_trdtrc )   ztrtrd(:,:,:) = ptrn(:,:,:,jn)                        ! save input trn for trend computation
269            !
270            WHERE( ptrn(:,:,:,jn) < 0. )   ptrn(:,:,:,jn) = 0.
271            !
272            IF( l_trdtrc ) THEN
273               ztrtrd(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt
274               CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radn, ztrtrd )       ! standard     trend handling
275            ENDIF
276            !
277         END DO
278         !
279      ENDIF
280      !
281      IF( l_trdtrc )  DEALLOCATE( ztrtrd )
282      !
283   END SUBROUTINE trc_rad_sms
284
285#else
286   !!----------------------------------------------------------------------
287   !!   Dummy module :                                         NO TOP model
288   !!----------------------------------------------------------------------
289CONTAINS
290   SUBROUTINE trc_rad( kt )              ! Empty routine
291      INTEGER, INTENT(in) ::   kt
292      WRITE(*,*) 'trc_rad: You should not have seen this print! error?', kt
293   END SUBROUTINE trc_rad
294#endif
295   
296   !!======================================================================
297END MODULE trcrad
Note: See TracBrowser for help on using the repository browser.