Changeset 10425 for NEMO/trunk/src/TOP/TRP/trcrad.F90
- Timestamp:
- 2018-12-19T22:54:16+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/TOP/TRP/trcrad.F90
r10068 r10425 19 19 USE trdtra 20 20 USE prtctl_trc ! Print control for debbuging 21 USE lib_fortran 21 22 22 23 IMPLICIT NONE … … 27 28 28 29 LOGICAL , PUBLIC :: ln_trcrad !: flag to artificially correct negative concentrations 30 REAL(wp), DIMENSION(:,:), ALLOCATABLE:: gainmass 29 31 30 32 !!---------------------------------------------------------------------- … … 104 106 ENDIF 105 107 ENDIF 108 ! 109 ALLOCATE( gainmass(jptra,2) ) 110 gainmass(:,:) = 0. 106 111 ! 107 112 END SUBROUTINE trc_rad_ini … … 129 134 CHARACTER( len = 1), OPTIONAL , INTENT(in ) :: cpreserv ! flag to preserve content or not 130 135 ! 131 INTEGER :: ji, jj, jk, jn ! dummy loop indices 132 LOGICAL :: lldebug = .FALSE. ! local logical 133 REAL(wp):: ztrcorb, ztrmasb, zs2rdt ! temporary scalars 134 REAL(wp):: zcoef , ztrcorn, ztrmasn ! - - 135 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrtrdb, ztrtrdn ! workspace arrays 136 !!---------------------------------------------------------------------- 137 ! 138 IF( l_trdtrc ) ALLOCATE( ztrtrdb(jpi,jpj,jpk), ztrtrdn(jpi,jpj,jpk) ) 136 INTEGER :: ji, ji2, jj, jj2, jk, jn ! dummy loop indices 137 INTEGER :: icnt 138 LOGICAL :: lldebug = .FALSE. ! local logical 139 REAL(wp):: zcoef, zs2rdt, ztotmass 140 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrneg, ztrpos 141 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrtrd ! workspace arrays 142 !!---------------------------------------------------------------------- 143 ! 144 IF( l_trdtrc ) ALLOCATE( ztrtrd(jpi,jpj,jpk) ) 145 zs2rdt = 1. / ( 2. * rdt * REAL( nn_dttrc, wp ) ) 139 146 ! 140 147 IF( PRESENT( cpreserv ) ) THEN !== total tracer concentration is preserved ==! 141 148 ! 142 DO jn = jp_sms0, jp_sms1 143 ! 144 ztrcorb = 0._wp ; ztrmasb = 0._wp 145 ztrcorn = 0._wp ; ztrmasn = 0._wp 146 ! 147 IF( l_trdtrc ) THEN 148 ztrtrdb(:,:,:) = ptrb(:,:,:,jn) ! save input trb for trend computation 149 ztrtrdn(:,:,:) = ptrn(:,:,:,jn) ! save input trn for trend computation 150 ENDIF 151 ! ! sum over the global domain 152 ztrcorb = glob_sum( MIN( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:) ) 153 ztrcorn = glob_sum( MIN( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:) ) 154 ! 155 ztrmasb = glob_sum( MAX( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:) ) 156 ztrmasn = glob_sum( MAX( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:) ) 157 ! 158 IF( ztrcorb /= 0 ) THEN 159 zcoef = 1. + ztrcorb / ztrmasb 160 DO jk = 1, jpkm1 161 ptrb(:,:,jk,jn) = MAX( 0., ptrb(:,:,jk,jn) ) 162 ptrb(:,:,jk,jn) = ptrb(:,:,jk,jn) * zcoef * tmask(:,:,jk) 163 END DO 164 ENDIF 165 ! 166 IF( ztrcorn /= 0 ) THEN 167 zcoef = 1. + ztrcorn / ztrmasn 168 DO jk = 1, jpkm1 169 ptrn(:,:,jk,jn) = MAX( 0., ptrn(:,:,jk,jn) ) 170 ptrn(:,:,jk,jn) = ptrn(:,:,jk,jn) * zcoef * tmask(:,:,jk) 171 END DO 172 ENDIF 173 ! 174 IF( l_trdtrc ) THEN 175 ! 176 zs2rdt = 1. / ( 2. * rdt ) 177 ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt 178 ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt 179 CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb ) ! Asselin-like trend handling 180 CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn ) ! standard trend handling 181 ! 182 ENDIF 183 ! 184 END DO 185 ! 186 ELSE !== total CFC content is NOT strictly preserved ==! 187 ! 188 DO jn = jp_sms0, jp_sms1 189 ! 190 IF( l_trdtrc ) THEN 191 ztrtrdb(:,:,:) = ptrb(:,:,:,jn) ! save input trb for trend computation 192 ztrtrdn(:,:,:) = ptrn(:,:,:,jn) ! save input trn for trend computation 193 ENDIF 149 ALLOCATE( ztrneg(1:jpi,1:jpj,jp_sms0:jp_sms1), ztrpos(1:jpi,1:jpj,jp_sms0:jp_sms1) ) 150 151 DO jn = jp_sms0, jp_sms1 152 ztrneg(:,:,jn) = SUM( MIN( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:), dim = 3 ) ! sum of the negative values 153 ztrpos(:,:,jn) = SUM( MAX( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:), dim = 3 ) ! sum of the positive values 154 END DO 155 CALL sum3x3( ztrneg ) 156 CALL sum3x3( ztrpos ) 157 158 DO jn = jp_sms0, jp_sms1 159 ! 160 IF( l_trdtrc ) ztrtrd(:,:,:) = ptrb(:,:,:,jn) ! save input trb for trend computation 194 161 ! 195 162 DO jk = 1, jpkm1 196 163 DO jj = 1, jpj 197 164 DO ji = 1, jpi 198 ptrn(ji,jj,jk,jn) = MAX( 0. , ptrn(ji,jj,jk,jn) ) 199 ptrb(ji,jj,jk,jn) = MAX( 0. , ptrb(ji,jj,jk,jn) ) 165 IF( ztrneg(ji,jj,jn) /= 0. ) THEN ! if negative values over the 3x3 box 166 ! 167 ptrb(ji,jj,jk,jn) = ptrb(ji,jj,jk,jn) * tmask(ji,jj,jk) ! really needed? 168 IF( ptrb(ji,jj,jk,jn) < 0. ) ptrb(ji,jj,jk,jn) = 0. ! supress negative values 169 IF( ptrb(ji,jj,jk,jn) > 0. ) THEN ! use positive values to compensate mass gain 170 zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn) ! ztrpos > 0 as ptrb > 0 171 ptrb(ji,jj,jk,jn) = ptrb(ji,jj,jk,jn) * zcoef 172 IF( zcoef < 0. ) THEN ! if the compensation exceed the positive value 173 gainmass(jn,1) = gainmass(jn,1) - ptrb(ji,jj,jk,jn) * cvol(ji,jj,jk) ! we are adding mass... 174 ptrb(ji,jj,jk,jn) = 0. ! limit the compensation to keep positive value 175 ENDIF 176 ENDIF 177 ! 178 ENDIF 200 179 END DO 201 180 END DO … … 203 182 ! 204 183 IF( l_trdtrc ) THEN 205 ! 206 zs2rdt = 1. / ( 2. * rdt * REAL( nn_dttrc, wp ) ) 207 ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt 208 ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt 209 CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb ) ! Asselin-like trend handling 210 CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn ) ! standard trend handling 211 ! 184 ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt 185 CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrd ) ! Asselin-like trend handling 186 ENDIF 187 ! 188 END DO 189 190 IF( kt == nitend ) THEN 191 CALL mpp_sum( 'trcrad', gainmass(:,1) ) 192 DO jn = jp_sms0, jp_sms1 193 IF( gainmass(jn,1) > 0. ) THEN 194 ztotmass = glob_sum( 'trcrad', ptrb(:,:,:,jn) * cvol(:,:,:) ) 195 IF(lwp) WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrb, traceur ', jn & 196 & , ' total mass : ', ztotmass, ', mass gain : ', gainmass(jn,1) 197 END IF 198 END DO 199 ENDIF 200 201 DO jn = jp_sms0, jp_sms1 202 ztrneg(:,:,jn) = SUM( MIN( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:), dim = 3 ) ! sum of the negative values 203 ztrpos(:,:,jn) = SUM( MAX( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:), dim = 3 ) ! sum of the positive values 204 END DO 205 CALL sum3x3( ztrneg ) 206 CALL sum3x3( ztrpos ) 207 208 DO jn = jp_sms0, jp_sms1 209 ! 210 IF( l_trdtrc ) ztrtrd(:,:,:) = ptrn(:,:,:,jn) ! save input trb for trend computation 211 ! 212 DO jk = 1, jpkm1 213 DO jj = 1, jpj 214 DO ji = 1, jpi 215 IF( ztrneg(ji,jj,jn) /= 0. ) THEN ! if negative values over the 3x3 box 216 ! 217 ptrn(ji,jj,jk,jn) = ptrn(ji,jj,jk,jn) * tmask(ji,jj,jk) ! really needed? 218 IF( ptrn(ji,jj,jk,jn) < 0. ) ptrn(ji,jj,jk,jn) = 0. ! supress negative values 219 IF( ptrn(ji,jj,jk,jn) > 0. ) THEN ! use positive values to compensate mass gain 220 zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn) ! ztrpos > 0 as ptrb > 0 221 ptrn(ji,jj,jk,jn) = ptrn(ji,jj,jk,jn) * zcoef 222 IF( zcoef < 0. ) THEN ! if the compensation exceed the positive value 223 gainmass(jn,2) = gainmass(jn,2) - ptrn(ji,jj,jk,jn) * cvol(ji,jj,jk) ! we are adding mass... 224 ptrn(ji,jj,jk,jn) = 0. ! limit the compensation to keep positive value 225 ENDIF 226 ENDIF 227 ! 228 ENDIF 229 END DO 230 END DO 231 END DO 232 ! 233 IF( l_trdtrc ) THEN 234 ztrtrd(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt 235 CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrd ) ! standard trend handling 236 ENDIF 237 ! 238 END DO 239 240 IF( kt == nitend ) THEN 241 CALL mpp_sum( 'trcrad', gainmass(:,2) ) 242 DO jn = jp_sms0, jp_sms1 243 IF( gainmass(jn,2) > 0. ) THEN 244 ztotmass = glob_sum( 'trcrad', ptrn(:,:,:,jn) * cvol(:,:,:) ) 245 WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrn, traceur ', jn & 246 & , ' total mass : ', ztotmass, ', mass gain : ', gainmass(jn,1) 247 END IF 248 END DO 249 ENDIF 250 251 DEALLOCATE( ztrneg, ztrpos ) 252 ! 253 ELSE !== total CFC content is NOT strictly preserved ==! 254 ! 255 DO jn = jp_sms0, jp_sms1 256 ! 257 IF( l_trdtrc ) ztrtrd(:,:,:) = ptrb(:,:,:,jn) ! save input trb for trend computation 258 ! 259 WHERE( ptrb(:,:,:,jn) < 0. ) ptrb(:,:,:,jn) = 0. 260 ! 261 IF( l_trdtrc ) THEN 262 ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt 263 CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrd ) ! Asselin-like trend handling 264 ENDIF 265 ! 266 IF( l_trdtrc ) ztrtrd(:,:,:) = ptrn(:,:,:,jn) ! save input trn for trend computation 267 ! 268 WHERE( ptrn(:,:,:,jn) < 0. ) ptrn(:,:,:,jn) = 0. 269 ! 270 IF( l_trdtrc ) THEN 271 ztrtrd(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt 272 CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrd ) ! standard trend handling 212 273 ENDIF 213 274 ! … … 216 277 ENDIF 217 278 ! 218 IF( l_trdtrc ) DEALLOCATE( ztrtrd b, ztrtrdn)279 IF( l_trdtrc ) DEALLOCATE( ztrtrd ) 219 280 ! 220 281 END SUBROUTINE trc_rad_sms
Note: See TracChangeset
for help on using the changeset viewer.