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 branches/UKMO/CO6_KD490/NEMOGCM/NEMO/TOP_SRC/TRP – NEMO

source: branches/UKMO/CO6_KD490/NEMOGCM/NEMO/TOP_SRC/TRP/trcrad.F90 @ 6332

Last change on this file since 6332 was 6332, checked in by deazer, 8 years ago

Tested Initial run one day physics only in rose suite.
File size: 11.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 oce_trc             ! ocean dynamics and tracers variables
16   USE trc                 ! ocean passive tracers variables
17   USE trd_oce
18   USE trdtra
19   USE prtctl_trc          ! Print control for debbuging
20#if defined key_tracer_budget
21   USE iom
22#endif
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC trc_rad         ! routine called by trcstp.F90
28
29   !! * Substitutions
30#  include "top_substitute.h90"
31   !!----------------------------------------------------------------------
32   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
33   !! $Id$
34   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
35   !!----------------------------------------------------------------------
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      CHARACTER (len=22) :: charout
56      !!----------------------------------------------------------------------
57      !
58      IF( nn_timing == 1 )  CALL timing_start('trc_rad')
59      !
60      IF( kt == nittrc000 ) THEN
61         IF(lwp) WRITE(numout,*)
62         IF(lwp) WRITE(numout,*) 'trc_rad : Correct artificial negative concentrations '
63         IF(lwp) WRITE(numout,*) '~~~~~~~ '
64      ENDIF
65
66      IF( lk_cfc     )   CALL trc_rad_sms( kt, trb, trn, jp_cfc0 , jp_cfc1               )  ! CFC model
67      IF( lk_c14b    )   CALL trc_rad_sms( kt, trb, trn, jp_c14b0, jp_c14b1              )  ! bomb C14
68      IF( lk_pisces  )   CALL trc_rad_sms( kt, trb, trn, jp_pcs0 , jp_pcs1, cpreserv='Y' )  ! PISCES model
69      IF( lk_my_trc  )   CALL trc_rad_sms( kt, trb, trn, jp_myt0 , jp_myt1               )  ! MY_TRC model
70
71      !
72      IF(ln_ctl) 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=trn, mask=tmask, clinfo=ctrcnm )
76      ENDIF
77      !
78      IF( nn_timing == 1 )  CALL timing_stop('trc_rad')
79      !
80   END SUBROUTINE trc_rad
81
82   SUBROUTINE trc_rad_sms( kt, ptrb, ptrn, jp_sms0, jp_sms1, cpreserv )
83      !!-----------------------------------------------------------------------------
84      !!                  ***  ROUTINE trc_rad_sms  ***
85      !!
86      !! ** Purpose :   "crappy" routine to correct artificial negative
87      !!              concentrations due to isopycnal scheme
88      !!
89      !! ** Method  : 2 cases :
90      !!                - Set negative concentrations to zero while computing
91      !!                  the corresponding tracer content that is added to the
92      !!                  tracers. Then, adjust the tracer concentration using
93      !!                  a multiplicative factor so that the total tracer
94      !!                  concentration is preserved.
95      !!                - simply set to zero the negative CFC concentration
96      !!                  (the total content of concentration is not strictly preserved)
97      !!--------------------------------------------------------------------------------
98      !! Arguments
99      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
100      INTEGER  , INTENT( in ) ::  &
101         jp_sms0, &       !: First index of the passive tracer model
102         jp_sms1          !: Last  index of  the passive tracer model
103
104      REAL(wp), DIMENSION (jpi,jpj,jpk,jptra), INTENT( inout )  :: &
105         ptrb, ptrn       !: before and now traceur concentration
106
107      CHARACTER( len = 1) , INTENT(in), OPTIONAL  :: &
108         cpreserv          !: flag to preserve content or not
109     
110      ! Local declarations
111      INTEGER  :: ji, jj, jk, jn     ! dummy loop indices
112      REAL(wp) :: ztrcorb, ztrmasb   ! temporary scalars
113      REAL(wp) :: zcoef, ztrcorn, ztrmasn   !    "         "
114      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrtrdb, ztrtrdn   ! workspace arrays
115#if defined key_tracer_budget
116      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  ztrtrdb_m1 ! slwa
117#endif
118      REAL(wp) :: zs2rdt
119      LOGICAL ::   lldebug = .FALSE.
120      !!----------------------------------------------------------------------
121
122 
123      IF( l_trdtrc )  CALL wrk_alloc( jpi, jpj, jpk, ztrtrdb, ztrtrdn )
124#if defined key_tracer_budget
125      IF( kt == nittrc000 .AND. l_trdtrc) THEN
126         ALLOCATE( ztrtrdb_m1(jpi,jpj,jpk,jptra) )  ! slwa
127         IF( ln_rsttr .AND.    &                     ! Restart: read in restart  file
128            iom_varid( numrtr, 'rdb_trend_'//TRIM(ctrcnm(1)), ldstop = .FALSE. ) > 0 ) THEN
129            IF(lwp) WRITE(numout,*) '          nittrc000-nn_dttrc RDB tracer trend read in the restart file'
130            DO jn = 1, jptra
131               CALL iom_get( numrtr, jpdom_autoglo, 'rdb_trend_'//TRIM(ctrcnm(jn)), ztrtrdb_m1(:,:,:,jn) )   ! before tracer trend for rdb
132            END DO
133         ELSE
134           ztrtrdb_m1=0.0
135         ENDIF
136      ENDIF
137#endif
138     
139      IF( PRESENT( cpreserv )  ) THEN   !  total tracer concentration is preserved
140     
141         DO jn = jp_sms0, jp_sms1
142            !                                                        ! ===========
143            ztrcorb = 0.e0   ;   ztrmasb = 0.e0
144            ztrcorn = 0.e0   ;   ztrmasn = 0.e0
145
146            IF( l_trdtrc ) THEN
147               ztrtrdb(:,:,:) = ptrb(:,:,:,jn)                        ! save input trb for trend computation
148               ztrtrdn(:,:,:) = ptrn(:,:,:,jn)                        ! save input trn for trend computation
149            ENDIF
150            !                                                         ! sum over the global domain
151            ztrcorb = glob_sum( MIN( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:) )
152            ztrcorn = glob_sum( MIN( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:) )
153
154            ztrmasb = glob_sum( MAX( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:) )
155            ztrmasn = glob_sum( MAX( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:) )
156
157            IF( ztrcorb /= 0 ) THEN
158               zcoef = 1. + ztrcorb / ztrmasb
159               DO jk = 1, jpkm1
160                  ptrb(:,:,jk,jn) = MAX( 0., ptrb(:,:,jk,jn) )
161                  ptrb(:,:,jk,jn) = ptrb(:,:,jk,jn) * zcoef * tmask(:,:,jk)
162               END DO
163            ENDIF
164
165            IF( ztrcorn /= 0 ) THEN
166               zcoef = 1. + ztrcorn / ztrmasn
167               DO jk = 1, jpkm1
168                  ptrn(:,:,jk,jn) = MAX( 0., ptrn(:,:,jk,jn) )
169                  ptrn(:,:,jk,jn) = ptrn(:,:,jk,jn) * zcoef * tmask(:,:,jk)
170               END DO
171            ENDIF
172            !
173            IF( l_trdtrc ) THEN
174               !
175               zs2rdt = 1. / ( 2. * rdt )
176               ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt
177               ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt 
178#if defined key_tracer_budget
179! slwa budget code
180               DO jk = 1, jpkm1
181                  ztrtrdb(:,:,jk) = ztrtrdb(:,:,jk) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk)
182                  ztrtrdn(:,:,jk) = ztrtrdn(:,:,jk) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk)
183               END DO
184               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb_m1(:,:,:,jn) )
185               ztrtrdb_m1(:,:,:,jn)=ztrtrdb(:,:,:)
186#else
187               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb )       ! Asselin-like trend handling
188#endif
189               CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn )       ! standard     trend handling
190              !
191            ENDIF
192
193         END DO
194         !
195         !
196      ELSE  ! total CFC content is not strictly preserved
197
198         DO jn = jp_sms0, jp_sms1 
199
200           IF( l_trdtrc ) THEN
201              ztrtrdb(:,:,:) = ptrb(:,:,:,jn)                        ! save input trb for trend computation
202              ztrtrdn(:,:,:) = ptrn(:,:,:,jn)                        ! save input trn for trend computation
203           ENDIF
204
205            DO jk = 1, jpkm1
206               DO jj = 1, jpj
207                  DO ji = 1, jpi
208                     ptrn(ji,jj,jk,jn) = MAX( 0. , ptrn(ji,jj,jk,jn) )
209                     ptrb(ji,jj,jk,jn) = MAX( 0. , ptrb(ji,jj,jk,jn) )
210                  END DO
211               END DO
212            END DO
213         
214            IF( l_trdtrc ) THEN
215               !
216               zs2rdt = 1. / ( 2. * rdt * FLOAT( nn_dttrc ) )
217               ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt
218               ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt 
219#if defined key_tracer_budget
220! slwa budget code
221               DO jk = 1, jpkm1
222                  ztrtrdb(:,:,jk) = ztrtrdb(:,:,jk) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk)
223                  ztrtrdn(:,:,jk) = ztrtrdn(:,:,jk) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk)
224               END DO
225               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb_m1(:,:,:,jn) )
226               ztrtrdb_m1(:,:,:,jn)=ztrtrdb(:,:,:)
227#else
228               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb )       ! Asselin-like trend handling
229#endif
230               CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn )       ! standard     trend handling
231              !
232            ENDIF
233            !
234         ENDDO
235
236      ENDIF
237
238#if defined key_tracer_budget
239      !                                           Write in the tracer restart file
240      !                                          *******************************
241      IF( lrst_trc ) THEN
242         IF(lwp) WRITE(numout,*)
243         IF(lwp) WRITE(numout,*) 'trc : RDB trend at last time step for tracer budget written in tracer restart file ',   &
244            &                    'at it= ', kt,' date= ', ndastp
245         IF(lwp) WRITE(numout,*) '~~~~'
246         DO jn = 1, jptra
247            CALL iom_rstput( kt, nitrst, numrtw, 'rdb_trend_'//TRIM(ctrcnm(jn)), ztrtrdb_m1(:,:,:,jn) )
248         END DO
249      ENDIF
250#endif
251
252      IF( l_trdtrc )  CALL wrk_dealloc( jpi, jpj, jpk, ztrtrdb, ztrtrdn )
253
254   END SUBROUTINE trc_rad_sms
255#else
256   !!----------------------------------------------------------------------
257   !!   Dummy module :                                         NO TOP model
258   !!----------------------------------------------------------------------
259CONTAINS
260   SUBROUTINE trc_rad( kt )              ! Empty routine
261      INTEGER, INTENT(in) ::   kt
262      WRITE(*,*) 'trc_rad: You should not have seen this print! error?', kt
263   END SUBROUTINE trc_rad
264#endif
265   
266   !!======================================================================
267END MODULE trcrad
Note: See TracBrowser for help on using the repository browser.