Changeset 11434 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src
- Timestamp:
- 2019-08-13T09:18:16+02:00 (5 years ago)
- Location:
- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP
- Files:
-
- 2 edited
-
trcrad.F90 (modified) (5 diffs)
-
trctrp.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trcrad.F90
r10985 r11434 37 37 CONTAINS 38 38 39 SUBROUTINE trc_rad( kt, Kbb, Kmm, Krhs,ptr )39 SUBROUTINE trc_rad( kt, Kbb, Kmm, ptr ) 40 40 !!---------------------------------------------------------------------- 41 41 !! *** ROUTINE trc_rad *** … … 52 52 !! (the total CFC content is not strictly preserved) 53 53 !!---------------------------------------------------------------------- 54 INTEGER, INTENT(in ) :: kt 55 INTEGER, INTENT(in ) :: Kbb, Kmm , Krhs! time level indices56 REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr 54 INTEGER, INTENT(in ) :: kt ! ocean time-step index 55 INTEGER, INTENT(in ) :: Kbb, Kmm ! time level indices 56 REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr ! passive tracers and RHS of tracer equation 57 57 ! 58 58 CHARACTER (len=22) :: charout … … 61 61 IF( ln_timing ) CALL timing_start('trc_rad') 62 62 ! 63 IF( ln_age ) CALL trc_rad_sms( kt, K mm, Krhs, ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Kmm), jp_age , jp_age ) ! AGE64 IF( ll_cfc ) CALL trc_rad_sms( kt, K mm, Krhs, ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Kmm), jp_cfc0, jp_cfc1 ) ! CFC model65 IF( ln_c14 ) CALL trc_rad_sms( kt, K mm, Krhs, ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Kmm), jp_c14 , jp_c14 ) ! C1466 IF( ln_pisces ) CALL trc_rad_sms( kt, K mm, Krhs, ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Kmm), jp_pcs0, jp_pcs1, cpreserv='Y' ) ! PISCES model67 IF( ln_my_trc ) CALL trc_rad_sms( kt, K mm, Krhs, ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Kmm), jp_myt0, jp_myt1 ) ! MY_TRC model63 IF( ln_age ) CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_age , jp_age ) ! AGE 64 IF( ll_cfc ) CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_cfc0, jp_cfc1 ) ! CFC model 65 IF( ln_c14 ) CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_c14 , jp_c14 ) ! C14 66 IF( ln_pisces ) CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_pcs0, jp_pcs1, cpreserv='Y' ) ! PISCES model 67 IF( ln_my_trc ) CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_myt0, jp_myt1 ) ! MY_TRC model 68 68 ! 69 69 IF(ln_ctl) THEN ! print mean trends (used for debugging) 70 70 WRITE(charout, FMT="('rad')") 71 71 CALL prt_ctl_trc_info( charout ) 72 CALL prt_ctl_trc( tab4d=ptr(:,:,:,:,K mm), mask=tmask, clinfo=ctrcnm )72 CALL prt_ctl_trc( tab4d=ptr(:,:,:,:,Kbb), mask=tmask, clinfo=ctrcnm ) 73 73 ENDIF 74 74 ! … … 115 115 116 116 117 SUBROUTINE trc_rad_sms( kt, Kmm, Krhs, ptrb, ptrn, jp_sms0, jp_sms1, cpreserv ) 118 !!----------------------------------------------------------------------------- 119 !! *** ROUTINE trc_rad_sms *** 120 !! 121 !! ** Purpose : "crappy" routine to correct artificial negative 122 !! concentrations due to isopycnal scheme 123 !! 124 !! ** Method : 2 cases : 125 !! - Set negative concentrations to zero while computing 126 !! the corresponding tracer content that is added to the 127 !! tracers. Then, adjust the tracer concentration using 128 !! a multiplicative factor so that the total tracer 129 !! concentration is preserved. 130 !! - simply set to zero the negative CFC concentration 131 !! (the total content of concentration is not strictly preserved) 132 !!-------------------------------------------------------------------------------- 133 INTEGER , INTENT(in ) :: kt ! ocean time-step index 134 INTEGER , INTENT(in ) :: Kmm, Krhs ! time level indices 135 INTEGER , INTENT(in ) :: jp_sms0, jp_sms1 ! First & last index of the passive tracer model 136 REAL(wp), DIMENSION (jpi,jpj,jpk,jptra), INTENT(inout) :: ptrb , ptrn ! before and now traceur concentration 137 CHARACTER( len = 1), OPTIONAL , INTENT(in ) :: cpreserv ! flag to preserve content or not 138 ! 139 INTEGER :: ji, ji2, jj, jj2, jk, jn ! dummy loop indices 140 INTEGER :: icnt 141 LOGICAL :: lldebug = .FALSE. ! local logical 142 REAL(wp):: zcoef, zs2rdt, ztotmass 143 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrneg, ztrpos 144 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrtrd ! workspace arrays 145 !!---------------------------------------------------------------------- 146 ! 147 IF( l_trdtrc ) ALLOCATE( ztrtrd(jpi,jpj,jpk) ) 148 zs2rdt = 1. / ( 2. * rdt * REAL( nn_dttrc, wp ) ) 149 ! 150 IF( PRESENT( cpreserv ) ) THEN !== total tracer concentration is preserved ==! 151 ! 152 ALLOCATE( ztrneg(1:jpi,1:jpj,jp_sms0:jp_sms1), ztrpos(1:jpi,1:jpj,jp_sms0:jp_sms1) ) 153 154 DO jn = jp_sms0, jp_sms1 155 ztrneg(:,:,jn) = SUM( MIN( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:), dim = 3 ) ! sum of the negative values 156 ztrpos(:,:,jn) = SUM( MAX( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:), dim = 3 ) ! sum of the positive values 157 END DO 158 CALL sum3x3( ztrneg ) 159 CALL sum3x3( ztrpos ) 160 161 DO jn = jp_sms0, jp_sms1 162 ! 163 IF( l_trdtrc ) ztrtrd(:,:,:) = ptrb(:,:,:,jn) ! save input tr(:,:,:,:,Kbb) for trend computation 164 ! 165 DO jk = 1, jpkm1 166 DO jj = 1, jpj 167 DO ji = 1, jpi 168 IF( ztrneg(ji,jj,jn) /= 0. ) THEN ! if negative values over the 3x3 box 169 ! 170 ptrb(ji,jj,jk,jn) = ptrb(ji,jj,jk,jn) * tmask(ji,jj,jk) ! really needed? 171 IF( ptrb(ji,jj,jk,jn) < 0. ) ptrb(ji,jj,jk,jn) = 0. ! supress negative values 172 IF( ptrb(ji,jj,jk,jn) > 0. ) THEN ! use positive values to compensate mass gain 173 zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn) ! ztrpos > 0 as ptrb > 0 174 ptrb(ji,jj,jk,jn) = ptrb(ji,jj,jk,jn) * zcoef 175 IF( zcoef < 0. ) THEN ! if the compensation exceed the positive value 176 gainmass(jn,1) = gainmass(jn,1) - ptrb(ji,jj,jk,jn) * cvol(ji,jj,jk) ! we are adding mass... 177 ptrb(ji,jj,jk,jn) = 0. ! limit the compensation to keep positive value 178 ENDIF 179 ENDIF 180 ! 181 ENDIF 182 END DO 183 END DO 184 END DO 185 ! 186 IF( l_trdtrc ) THEN 187 ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt 188 CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radb, ztrtrd ) ! Asselin-like trend handling 189 ENDIF 190 ! 191 END DO 192 193 IF( kt == nitend ) THEN 194 CALL mpp_sum( 'trcrad', gainmass(:,1) ) 195 DO jn = jp_sms0, jp_sms1 196 IF( gainmass(jn,1) > 0. ) THEN 197 ztotmass = glob_sum( 'trcrad', ptrb(:,:,:,jn) * cvol(:,:,:) ) 198 IF(lwp) WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrb, traceur ', jn & 199 & , ' total mass : ', ztotmass, ', mass gain : ', gainmass(jn,1) 200 END IF 201 END DO 202 ENDIF 203 204 DO jn = jp_sms0, jp_sms1 205 ztrneg(:,:,jn) = SUM( MIN( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:), dim = 3 ) ! sum of the negative values 206 ztrpos(:,:,jn) = SUM( MAX( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:), dim = 3 ) ! sum of the positive values 207 END DO 208 CALL sum3x3( ztrneg ) 209 CALL sum3x3( ztrpos ) 210 211 DO jn = jp_sms0, jp_sms1 212 ! 213 IF( l_trdtrc ) ztrtrd(:,:,:) = ptrn(:,:,:,jn) ! save input tr for trend computation 214 ! 215 DO jk = 1, jpkm1 216 DO jj = 1, jpj 217 DO ji = 1, jpi 218 IF( ztrneg(ji,jj,jn) /= 0. ) THEN ! if negative values over the 3x3 box 219 ! 220 ptrn(ji,jj,jk,jn) = ptrn(ji,jj,jk,jn) * tmask(ji,jj,jk) ! really needed? 221 IF( ptrn(ji,jj,jk,jn) < 0. ) ptrn(ji,jj,jk,jn) = 0. ! supress negative values 222 IF( ptrn(ji,jj,jk,jn) > 0. ) THEN ! use positive values to compensate mass gain 223 zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn) ! ztrpos > 0 as ptrb > 0 224 ptrn(ji,jj,jk,jn) = ptrn(ji,jj,jk,jn) * zcoef 225 IF( zcoef < 0. ) THEN ! if the compensation exceed the positive value 226 gainmass(jn,2) = gainmass(jn,2) - ptrn(ji,jj,jk,jn) * cvol(ji,jj,jk) ! we are adding mass... 227 ptrn(ji,jj,jk,jn) = 0. ! limit the compensation to keep positive value 228 ENDIF 229 ENDIF 230 ! 231 ENDIF 232 END DO 233 END DO 234 END DO 235 ! 236 IF( l_trdtrc ) THEN 237 ztrtrd(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt 238 CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radn, ztrtrd ) ! standard trend handling 239 ENDIF 240 ! 241 END DO 242 243 IF( kt == nitend ) THEN 244 CALL mpp_sum( 'trcrad', gainmass(:,2) ) 245 DO jn = jp_sms0, jp_sms1 246 IF( gainmass(jn,2) > 0. ) THEN 247 ztotmass = glob_sum( 'trcrad', ptrn(:,:,:,jn) * cvol(:,:,:) ) 248 WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrn, traceur ', jn & 249 & , ' total mass : ', ztotmass, ', mass gain : ', gainmass(jn,1) 250 END IF 251 END DO 252 ENDIF 253 254 DEALLOCATE( ztrneg, ztrpos ) 255 ! 256 ELSE !== total CFC content is NOT strictly preserved ==! 257 ! 258 DO jn = jp_sms0, jp_sms1 259 ! 260 IF( l_trdtrc ) ztrtrd(:,:,:) = ptrb(:,:,:,jn) ! save input tr for trend computation 261 ! 262 WHERE( ptrb(:,:,:,jn) < 0. ) ptrb(:,:,:,jn) = 0. 263 ! 264 IF( l_trdtrc ) THEN 265 ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt 266 CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radb, ztrtrd ) ! Asselin-like trend handling 267 ENDIF 268 ! 269 IF( l_trdtrc ) ztrtrd(:,:,:) = ptrn(:,:,:,jn) ! save input tr for trend computation 270 ! 271 WHERE( ptrn(:,:,:,jn) < 0. ) ptrn(:,:,:,jn) = 0. 272 ! 273 IF( l_trdtrc ) THEN 274 ztrtrd(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt 275 CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radn, ztrtrd ) ! standard trend handling 276 ENDIF 277 ! 278 END DO 279 ! 280 ENDIF 117 SUBROUTINE trc_rad_sms( kt, Kbb, Kmm, ptr, jp_sms0, jp_sms1, cpreserv ) 118 !!----------------------------------------------------------------------------- 119 !! *** ROUTINE trc_rad_sms *** 120 !! 121 !! ** Purpose : "crappy" routine to correct artificial negative 122 !! concentrations due to isopycnal scheme 123 !! 124 !! ** Method : 2 cases : 125 !! - Set negative concentrations to zero while computing 126 !! the corresponding tracer content that is added to the 127 !! tracers. Then, adjust the tracer concentration using 128 !! a multiplicative factor so that the total tracer 129 !! concentration is preserved. 130 !! - simply set to zero the negative CFC concentration 131 !! (the total content of concentration is not strictly preserved) 132 !!-------------------------------------------------------------------------------- 133 INTEGER , INTENT(in ) :: kt ! ocean time-step index 134 INTEGER , INTENT(in ) :: Kbb, Kmm ! time level indices 135 INTEGER , INTENT(in ) :: jp_sms0, jp_sms1 ! First & last index of the passive tracer model 136 REAL(wp), DIMENSION (jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr ! before and now traceur concentration 137 CHARACTER( len = 1), OPTIONAL , INTENT(in ) :: cpreserv ! flag to preserve content or not 138 ! 139 INTEGER :: ji, ji2, jj, jj2, jk, jn, jt ! dummy loop indices 140 INTEGER :: icnt, itime 141 LOGICAL :: lldebug = .FALSE. ! local logical 142 REAL(wp):: zcoef, zs2rdt, ztotmass 143 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrneg, ztrpos 144 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrtrd ! workspace arrays 145 !!---------------------------------------------------------------------- 146 ! 147 IF( l_trdtrc ) ALLOCATE( ztrtrd(jpi,jpj,jpk) ) 148 zs2rdt = 1. / ( 2. * rdt * REAL( nn_dttrc, wp ) ) 149 ! 150 DO jt = 1,2 ! Loop over time indices since exactly the same fix is applied to "now" and "after" fields 151 IF( jt == 1 ) itime = Kbb 152 IF( jt == 2 ) itime = Kmm 153 154 IF( PRESENT( cpreserv ) ) THEN !== total tracer concentration is preserved ==! 155 ! 156 ALLOCATE( ztrneg(1:jpi,1:jpj,jp_sms0:jp_sms1), ztrpos(1:jpi,1:jpj,jp_sms0:jp_sms1) ) 157 158 DO jn = jp_sms0, jp_sms1 159 ztrneg(:,:,jn) = SUM( MIN( 0., ptr(:,:,:,jn,itime) ) * cvol(:,:,:), dim = 3 ) ! sum of the negative values 160 ztrpos(:,:,jn) = SUM( MAX( 0., ptr(:,:,:,jn,itime) ) * cvol(:,:,:), dim = 3 ) ! sum of the positive values 161 END DO 162 CALL sum3x3( ztrneg ) 163 CALL sum3x3( ztrpos ) 164 165 DO jn = jp_sms0, jp_sms1 166 ! 167 IF( l_trdtrc ) ztrtrd(:,:,:) = ptr(:,:,:,jn,itime) ! save input tr(:,:,:,:,Kbb) for trend computation 168 ! 169 DO jk = 1, jpkm1 170 DO jj = 1, jpj 171 DO ji = 1, jpi 172 IF( ztrneg(ji,jj,jn) /= 0. ) THEN ! if negative values over the 3x3 box 173 ! 174 ptr(ji,jj,jk,jn,itime) = ptr(ji,jj,jk,jn,itime) * tmask(ji,jj,jk) ! really needed? 175 IF( ptr(ji,jj,jk,jn,itime) < 0. ) ptr(ji,jj,jk,jn,itime) = 0. ! suppress negative values 176 IF( ptr(ji,jj,jk,jn,itime) > 0. ) THEN ! use positive values to compensate mass gain 177 zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn) ! ztrpos > 0 as ptr > 0 178 ptr(ji,jj,jk,jn,itime) = ptr(ji,jj,jk,jn,itime) * zcoef 179 IF( zcoef < 0. ) THEN ! if the compensation exceed the positive value 180 gainmass(jn,1) = gainmass(jn,1) - ptr(ji,jj,jk,jn,itime) * cvol(ji,jj,jk) ! we are adding mass... 181 ptr(ji,jj,jk,jn,itime) = 0. ! limit the compensation to keep positive value 182 ENDIF 183 ENDIF 184 ! 185 ENDIF 186 END DO 187 END DO 188 END DO 189 ! 190 IF( l_trdtrc ) THEN 191 ztrtrd(:,:,:) = ( ptr(:,:,:,jn,itime) - ztrtrd(:,:,:) ) * zs2rdt 192 CALL trd_tra( kt, Kbb, Kmm, 'TRC', jn, jptra_radb, ztrtrd ) ! Asselin-like trend handling 193 ENDIF 194 ! 195 END DO 196 197 IF( kt == nitend ) THEN 198 CALL mpp_sum( 'trcrad', gainmass(:,1) ) 199 DO jn = jp_sms0, jp_sms1 200 IF( gainmass(jn,1) > 0. ) THEN 201 ztotmass = glob_sum( 'trcrad', ptr(:,:,:,jn,itime) * cvol(:,:,:) ) 202 IF(lwp) WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrb, traceur ', jn & 203 & , ' total mass : ', ztotmass, ', mass gain : ', gainmass(jn,1) 204 END IF 205 END DO 206 ENDIF 207 208 DEALLOCATE( ztrneg, ztrpos ) 209 ! 210 ELSE !== total CFC content is NOT strictly preserved ==! 211 ! 212 DO jn = jp_sms0, jp_sms1 213 ! 214 IF( l_trdtrc ) ztrtrd(:,:,:) = ptr(:,:,:,jn,itime) ! save input tr for trend computation 215 ! 216 WHERE( ptr(:,:,:,jn,itime) < 0. ) ptr(:,:,:,jn,itime) = 0. 217 ! 218 IF( l_trdtrc ) THEN 219 ztrtrd(:,:,:) = ( ptr(:,:,:,jn,itime) - ztrtrd(:,:,:) ) * zs2rdt 220 CALL trd_tra( kt, Kbb, Kmm, 'TRC', jn, jptra_radb, ztrtrd ) ! Asselin-like trend handling 221 ENDIF 222 ! 223 END DO 224 ! 225 ENDIF 226 ! 227 END DO 281 228 ! 282 229 IF( l_trdtrc ) DEALLOCATE( ztrtrd ) … … 289 236 !!---------------------------------------------------------------------- 290 237 CONTAINS 291 SUBROUTINE trc_rad( kt, Kbb, Kmm , Krhs) ! Empty routine238 SUBROUTINE trc_rad( kt, Kbb, Kmm ) ! Empty routine 292 239 INTEGER, INTENT(in) :: kt 293 INTEGER, INTENT(in) :: Kbb, Kmm , Krhs! time level indices240 INTEGER, INTENT(in) :: Kbb, Kmm ! time level indices 294 241 WRITE(*,*) 'trc_rad: You should not have seen this print! error?', kt 295 242 END SUBROUTINE trc_rad -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trctrp.F90
r11427 r11434 36 36 PUBLIC trc_trp ! called by trc_stp 37 37 38 INTEGER, SAVE :: N_save ! Save value of time index for time swapping for ln_top_euler=.true.39 40 38 !!---------------------------------------------------------------------- 41 39 !! NEMO/TOP 4.0 , NEMO Consortium (2018) … … 61 59 IF( ln_timing ) CALL timing_start('trc_trp') 62 60 ! 63 IF ( kt == nit000 ) N_save = Kbb64 !65 61 IF( .NOT. lk_c1d ) THEN 66 62 ! … … 85 81 ! 86 82 ! Swap TOP time levels (= Nrhs_trc, Nbb_trc etc) 87 IF( ln_top_euler ) THEN 88 ! For Euler timestepping we need the "before" and "now" fields to be the same. 89 ! Use N_save to ensure that the indices stay in sync with the (leapfrogging) OCE time indices for nn_dttrc=1. 90 Krhs = N_save 91 N_save = Kmm 92 Kbb = Kaa 93 ELSE 94 Krhs = Kbb 95 Kbb = Kmm 96 ENDIF 83 Krhs = Kbb 84 Kbb = Kmm 97 85 Kmm = Kaa 98 86 Kaa = Krhs 99 87 ! 100 IF( ln_trcrad ) CALL trc_rad ( kt, Kbb, Kmm, Krhs, tr) ! Correct artificial negative concentrations101 IF( ln_trcdmp_clo ) CALL trc_dmp_clo( kt, Kbb, Kmm ) 88 IF( ln_trcrad ) CALL trc_rad ( kt, Kbb, Kmm, tr ) ! Correct artificial negative concentrations 89 IF( ln_trcdmp_clo ) CALL trc_dmp_clo( kt, Kbb, Kmm ) ! internal damping trends on closed seas only 102 90 103 91 ! … … 114 102 Kaa = Krhs 115 103 ! 116 IF( ln_trcrad ) CALL trc_rad( kt, Kbb, Kmm, Krhs,tr ) ! Correct artificial negative concentrations104 IF( ln_trcrad ) CALL trc_rad( kt, Kbb, Kmm, tr ) ! Correct artificial negative concentrations 117 105 ! 118 106 END IF
Note: See TracChangeset
for help on using the changeset viewer.
