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/trunk/src/TOP/TRP – NEMO

source: NEMO/trunk/src/TOP/TRP/trcrad.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

  • Property svn:keywords set to Id
File size: 12.1 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   !!            4.1  !  08-19  (A. Coward, D. Storkey) tidy up using new time-level indices
9   !!----------------------------------------------------------------------
10#if defined key_top
11   !!----------------------------------------------------------------------
12   !!   'key_top'                                                TOP models
13   !!----------------------------------------------------------------------
14   !!   trc_rad    : correction of negative concentrations
15   !!----------------------------------------------------------------------
16   USE par_trc             ! need jptra, number of passive tracers
17   USE oce_trc             ! ocean dynamics and tracers variables
18   USE trc                 ! ocean passive tracers variables
19   USE trd_oce
20   USE trdtra
21   USE prtctl_trc          ! Print control for debbuging
22   USE lib_fortran
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC trc_rad     
28   PUBLIC trc_rad_ini 
29
30   LOGICAL , PUBLIC ::   ln_trcrad           !: flag to artificially correct negative concentrations
31   REAL(wp), DIMENSION(:,:), ALLOCATABLE::   gainmass
32
33   !! * Substitutions
34#  include "do_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
37   !! $Id$
38   !! Software governed by the CeCILL license (see ./LICENSE)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE trc_rad( kt, Kbb, Kmm, ptr )
43      !!----------------------------------------------------------------------
44      !!                  ***  ROUTINE trc_rad  ***
45      !!
46      !! ** Purpose :   "crappy" routine to correct artificial negative
47      !!              concentrations due to isopycnal scheme
48      !!
49      !! ** Method  : - PISCES or LOBSTER: Set negative concentrations to zero
50      !!                while computing the corresponding tracer content that
51      !!                is added to the tracers. Then, adjust the tracer
52      !!                concentration using a multiplicative factor so that
53      !!                the total tracer concentration is preserved.
54      !!              - CFC: simply set to zero the negative CFC concentration
55      !!                (the total CFC content is not strictly preserved)
56      !!----------------------------------------------------------------------
57      INTEGER,                                    INTENT(in   ) :: kt         ! ocean time-step index
58      INTEGER,                                    INTENT(in   ) :: Kbb, Kmm   ! time level indices
59      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr        ! passive tracers and RHS of tracer equation
60      !
61      CHARACTER (len=22) :: charout
62      !!----------------------------------------------------------------------
63      !
64      IF( ln_timing )   CALL timing_start('trc_rad')
65      !
66      IF( ln_age     )   CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_age , jp_age                )  !  AGE
67      IF( ll_cfc     )   CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_cfc0, jp_cfc1               )  !  CFC model
68      IF( ln_c14     )   CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_c14 , jp_c14                )  !  C14
69      IF( ln_pisces  )   CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_pcs0, jp_pcs1, cpreserv='Y' )  !  PISCES model
70      IF( ln_my_trc  )   CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_myt0, jp_myt1               )  !  MY_TRC model
71      !
72      IF(sn_cfctl%l_prttrc) THEN      ! print mean trends (used for debugging)
73         WRITE(charout, FMT="('rad')")
74         CALL prt_ctl_trc_info( charout )
75         CALL prt_ctl_trc( tab4d=ptr(:,:,:,:,Kbb), mask=tmask, clinfo=ctrcnm )
76      ENDIF
77      !
78      IF( ln_timing )   CALL timing_stop('trc_rad')
79      !
80   END SUBROUTINE trc_rad
81
82
83   SUBROUTINE trc_rad_ini
84      !!---------------------------------------------------------------------
85      !!                  ***  ROUTINE trc _rad_ini ***
86      !!
87      !! ** Purpose :   read  namelist options
88      !!----------------------------------------------------------------------
89      INTEGER ::   ios   ! Local integer output status for namelist read
90      !!
91      NAMELIST/namtrc_rad/ ln_trcrad
92      !!----------------------------------------------------------------------
93      !
94      READ  ( numnat_ref, namtrc_rad, IOSTAT = ios, ERR = 907)
95907   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtrc_rad in reference namelist' )
96      READ  ( numnat_cfg, namtrc_rad, IOSTAT = ios, ERR = 908 )
97908   IF( ios > 0 )   CALL ctl_nam ( ios , 'namtrc_rad in configuration namelist' )
98      IF(lwm) WRITE( numont, namtrc_rad )
99
100      IF(lwp) THEN                     !   ! Control print
101         WRITE(numout,*)
102         WRITE(numout,*) 'trc_rad : Correct artificial negative concentrations '
103         WRITE(numout,*) '~~~~~~~ '
104         WRITE(numout,*) '   Namelist namtrc_rad : treatment of negative concentrations'
105         WRITE(numout,*) '      correct artificially negative concen. or not   ln_trcrad = ', ln_trcrad
106         WRITE(numout,*)
107         IF( ln_trcrad ) THEN   ;   WRITE(numout,*) '      ===>>   ensure the global tracer conservation'
108         ELSE                   ;   WRITE(numout,*) '      ===>>   NO strict global tracer conservation'     
109         ENDIF
110      ENDIF
111      !
112      ALLOCATE( gainmass(jptra,2) )
113      gainmass(:,:) = 0.
114      !
115   END SUBROUTINE trc_rad_ini
116
117
118   SUBROUTINE trc_rad_sms( kt, Kbb, Kmm, ptr, jp_sms0, jp_sms1, cpreserv )
119     !!-----------------------------------------------------------------------------
120     !!                  ***  ROUTINE trc_rad_sms  ***
121     !!
122     !! ** Purpose :   "crappy" routine to correct artificial negative
123     !!              concentrations due to isopycnal scheme
124     !!
125     !! ** Method  : 2 cases :
126     !!                - Set negative concentrations to zero while computing
127     !!                  the corresponding tracer content that is added to the
128     !!                  tracers. Then, adjust the tracer concentration using
129     !!                  a multiplicative factor so that the total tracer
130     !!                  concentration is preserved.
131     !!                - simply set to zero the negative CFC concentration
132     !!                  (the total content of concentration is not strictly preserved)
133     !!--------------------------------------------------------------------------------
134     INTEGER                                    , INTENT(in   ) ::   kt                 ! ocean time-step index
135     INTEGER                                    , INTENT(in   ) ::   Kbb, Kmm           ! time level indices
136     INTEGER                                    , INTENT(in   ) ::   jp_sms0, jp_sms1   ! First & last index of the passive tracer model
137     REAL(wp), DIMENSION (jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::   ptr                ! before and now traceur concentration
138     CHARACTER( len = 1), OPTIONAL              , INTENT(in   ) ::   cpreserv           ! flag to preserve content or not
139     !
140     INTEGER ::   ji, ji2, jj, jj2, jk, jn, jt ! dummy loop indices
141     INTEGER ::   icnt, itime
142     LOGICAL ::   lldebug = .FALSE.            ! local logical
143     REAL(wp)::   zcoef, zs2rdt, ztotmass
144     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrneg, ztrpos
145     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrtrd   ! workspace arrays
146     !!----------------------------------------------------------------------
147     !
148     IF( l_trdtrc )   ALLOCATE( ztrtrd(jpi,jpj,jpk) )
149     zs2rdt = 1. / ( 2. * rdt )
150     !
151     DO jt = 1,2  ! Loop over time indices since exactly the same fix is applied to "now" and "after" fields
152        IF( jt == 1 ) itime = Kbb
153        IF( jt == 2 ) itime = Kmm
154
155        IF( PRESENT( cpreserv )  ) THEN     !==  total tracer concentration is preserved  ==!
156           !
157           ALLOCATE( ztrneg(1:jpi,1:jpj,jp_sms0:jp_sms1), ztrpos(1:jpi,1:jpj,jp_sms0:jp_sms1) )
158
159           DO jn = jp_sms0, jp_sms1
160              ztrneg(:,:,jn) = SUM( MIN( 0., ptr(:,:,:,jn,itime) ) * cvol(:,:,:), dim = 3 )   ! sum of the negative values
161              ztrpos(:,:,jn) = SUM( MAX( 0., ptr(:,:,:,jn,itime) ) * cvol(:,:,:), dim = 3 )   ! sum of the positive values
162           END DO
163           CALL sum3x3( ztrneg )
164           CALL sum3x3( ztrpos )
165
166           DO jn = jp_sms0, jp_sms1
167              !
168              IF( l_trdtrc )   ztrtrd(:,:,:) = ptr(:,:,:,jn,itime)                       ! save input tr(:,:,:,:,Kbb) for trend computation           
169              !
170              DO_3D_11_11( 1, jpkm1 )
171                 IF( ztrneg(ji,jj,jn) /= 0. ) THEN                                 ! if negative values over the 3x3 box
172                    !
173                    ptr(ji,jj,jk,jn,itime) = ptr(ji,jj,jk,jn,itime) * tmask(ji,jj,jk)   ! really needed?
174                    IF( ptr(ji,jj,jk,jn,itime) < 0. ) ptr(ji,jj,jk,jn,itime) = 0.       ! suppress negative values
175                    IF( ptr(ji,jj,jk,jn,itime) > 0. ) THEN                    ! use positive values to compensate mass gain
176                       zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn)       ! ztrpos > 0 as ptr > 0
177                       ptr(ji,jj,jk,jn,itime) = ptr(ji,jj,jk,jn,itime) * zcoef
178                       IF( zcoef < 0. ) THEN                                  ! if the compensation exceed the positive value
179                          gainmass(jn,1) = gainmass(jn,1) - ptr(ji,jj,jk,jn,itime) * cvol(ji,jj,jk)   ! we are adding mass...
180                          ptr(ji,jj,jk,jn,itime) = 0.                         ! limit the compensation to keep positive value
181                       ENDIF
182                    ENDIF
183                    !
184                 ENDIF
185              END_3D
186              !
187              IF( l_trdtrc ) THEN
188                 ztrtrd(:,:,:) = ( ptr(:,:,:,jn,itime) - ztrtrd(:,:,:) ) * zs2rdt
189                 CALL trd_tra( kt, Kbb, Kmm, 'TRC', jn, jptra_radb, ztrtrd )       ! Asselin-like trend handling
190              ENDIF
191              !
192           END DO
193
194           IF( kt == nitend ) THEN
195              CALL mpp_sum( 'trcrad', gainmass(:,1) )
196              DO jn = jp_sms0, jp_sms1
197                 IF( gainmass(jn,1) > 0. ) THEN
198                    ztotmass = glob_sum( 'trcrad', ptr(:,:,:,jn,itime) * cvol(:,:,:) )
199                    IF(lwp) WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrb, traceur ', jn  &
200                         &        , ' total mass : ', ztotmass, ', mass gain : ',  gainmass(jn,1)
201                 END IF
202              END DO
203           ENDIF
204
205           DEALLOCATE( ztrneg, ztrpos )
206           !
207        ELSE                                !==  total CFC content is NOT strictly preserved  ==!
208           !
209           DO jn = jp_sms0, jp_sms1 
210              !
211              IF( l_trdtrc )   ztrtrd(:,:,:) = ptr(:,:,:,jn,itime)                 ! save input tr for trend computation
212              !
213              WHERE( ptr(:,:,:,jn,itime) < 0. )   ptr(:,:,:,jn,itime) = 0.
214              !
215              IF( l_trdtrc ) THEN
216                 ztrtrd(:,:,:) = ( ptr(:,:,:,jn,itime) - ztrtrd(:,:,:) ) * zs2rdt
217                 CALL trd_tra( kt, Kbb, Kmm, 'TRC', jn, jptra_radb, ztrtrd )       ! Asselin-like trend handling
218              ENDIF
219              !
220           END DO
221           !
222        ENDIF
223        !
224      END DO
225      !
226      IF( l_trdtrc )  DEALLOCATE( ztrtrd )
227      !
228   END SUBROUTINE trc_rad_sms
229
230#else
231   !!----------------------------------------------------------------------
232   !!   Dummy module :                                         NO TOP model
233   !!----------------------------------------------------------------------
234CONTAINS
235   SUBROUTINE trc_rad( kt, Kbb, Kmm )              ! Empty routine
236      INTEGER, INTENT(in) ::   kt
237      INTEGER, INTENT(in) ::   Kbb, Kmm  ! time level indices
238      WRITE(*,*) 'trc_rad: You should not have seen this print! error?', kt
239   END SUBROUTINE trc_rad
240#endif
241   
242   !!======================================================================
243END MODULE trcrad
Note: See TracBrowser for help on using the repository browser.