Changeset 10408
- Timestamp:
- 2018-12-18T12:21:21+01:00 (6 years ago)
- Location:
- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src
- Files:
-
- 2 edited
-
OCE/lib_fortran.F90 (modified) (4 diffs)
-
TOP/TRP/trcrad.F90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/lib_fortran.F90
r10397 r10408 21 21 USE in_out_manager ! I/O manager 22 22 USE lib_mpp ! distributed memory computing 23 USE lbclnk ! ocean lateral boundary conditions 23 24 24 25 IMPLICIT NONE … … 28 29 PUBLIC glob_sum_full ! used in many places (masked with tmask_h, ie only over the halos) 29 30 PUBLIC local_sum ! used in trcrad, local operation before glob_sum_delay 31 PUBLIC sum3x3 ! used in trcrad, do a sum over 3x3 boxes 30 32 PUBLIC DDPDD ! also used in closea module 31 33 PUBLIC glob_min, glob_max … … 42 44 INTERFACE local_sum 43 45 MODULE PROCEDURE local_sum_2d, local_sum_3d 46 END INTERFACE 47 INTERFACE sum3x3 48 MODULE PROCEDURE sum3x3_2d, sum3x3_3d 44 49 END INTERFACE 45 50 INTERFACE glob_min … … 194 199 END FUNCTION local_sum_3d 195 200 201 ! ! FUNCTION sum3x3 ! 202 203 SUBROUTINE sum3x3_2d( p2d ) 204 !!----------------------------------------------------------------------- 205 !! *** routine sum3x3_2d *** 206 !! 207 !! ** Purpose : sum over 3x3 boxes 208 !!---------------------------------------------------------------------- 209 REAL(wp), DIMENSION (:,:), INTENT(inout) :: p2d 210 ! 211 INTEGER :: ji, ji2, jj, jj2 ! dummy loop indices 212 !!---------------------------------------------------------------------- 213 ! 214 IF( SIZE(p2d,1) /= jpi ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_2d, the first dimension is not equal to jpi' ) 215 IF( SIZE(p2d,2) /= jpj ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_2d, the second dimension is not equal to jpj' ) 216 ! 217 DO jj = 1, jpj 218 DO ji = 1, jpi 219 IF( MOD(mig(ji), 3) == 1 .AND. MOD(mjg(jj), 3) == 1 ) THEN ! bottom left corber of a 3x3 box 220 ji2 = MIN(mig(ji)+2, jpiglo) - nimpp + 1 ! right position of the box 221 jj2 = MIN(mjg(jj)+2, jpjglo) - njmpp + 1 ! upper position of the box 222 IF( ji2 <= jpi .AND. jj2 <= jpj ) THEN ! the box is fully included in the local mpi domain 223 p2d(ji:ji2,jj:jj2) = SUM(p2d(ji:ji2,jj:jj2)) 224 ENDIF 225 ENDIF 226 END DO 227 END DO 228 CALL lbc_lnk( 'lib_fortran', p2d, 'T', 1. ) 229 IF( nbondi /= -1 ) THEN 230 IF( MOD(mig( 1), 3) == 1 ) p2d( 1,:) = p2d( 2,:) 231 IF( MOD(mig( 1), 3) == 2 ) p2d( 2,:) = p2d( 1,:) 232 ENDIF 233 IF( nbondi /= 1 ) THEN 234 IF( MOD(mig(jpi-2), 3) == 1 ) p2d( jpi,:) = p2d(jpi-1,:) 235 IF( MOD(mig(jpi-2), 3) == 0 ) p2d(jpi-1,:) = p2d( jpi,:) 236 ENDIF 237 IF( nbondj /= -1 ) THEN 238 IF( MOD(mjg( 1), 3) == 1 ) p2d(:, 1) = p2d(:, 2) 239 IF( MOD(mjg( 1), 3) == 2 ) p2d(:, 2) = p2d(:, 1) 240 ENDIF 241 IF( nbondj /= 1 ) THEN 242 IF( MOD(mjg(jpj-2), 3) == 1 ) p2d(:, jpj) = p2d(:,jpj-1) 243 IF( MOD(mjg(jpj-2), 3) == 0 ) p2d(:,jpj-1) = p2d(:, jpj) 244 ENDIF 245 CALL lbc_lnk( 'lib_fortran', p2d, 'T', 1. ) 246 247 END SUBROUTINE sum3x3_2d 248 249 SUBROUTINE sum3x3_3d( p3d ) 250 !!----------------------------------------------------------------------- 251 !! *** routine sum3x3_3d *** 252 !! 253 !! ** Purpose : sum over 3x3 boxes 254 !!---------------------------------------------------------------------- 255 REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: p3d 256 ! 257 INTEGER :: ji, ji2, jj, jj2, jn ! dummy loop indices 258 INTEGER :: ipn ! Third dimension size 259 !!---------------------------------------------------------------------- 260 ! 261 IF( SIZE(p3d,1) /= jpi ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_3d, the first dimension is not equal to jpi' ) 262 IF( SIZE(p3d,2) /= jpj ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_3d, the second dimension is not equal to jpj' ) 263 ipn = SIZE(p3d,3) 264 ! 265 DO jn = 1, ipn 266 DO jj = 1, jpj 267 DO ji = 1, jpi 268 IF( MOD(mig(ji), 3) == 1 .AND. MOD(mjg(jj), 3) == 1 ) THEN ! bottom left corber of a 3x3 box 269 ji2 = MIN(mig(ji)+2, jpiglo) - nimpp + 1 ! right position of the box 270 jj2 = MIN(mjg(jj)+2, jpjglo) - njmpp + 1 ! upper position of the box 271 IF( ji2 <= jpi .AND. jj2 <= jpj ) THEN ! the box is fully included in the local mpi domain 272 p3d(ji:ji2,jj:jj2,jn) = SUM(p3d(ji:ji2,jj:jj2,jn)) 273 ENDIF 274 ENDIF 275 END DO 276 END DO 277 END DO 278 CALL lbc_lnk( 'lib_fortran', p3d, 'T', 1. ) 279 IF( nbondi /= -1 ) THEN 280 IF( MOD(mig( 1), 3) == 1 ) p3d( 1,:,:) = p3d( 2,:,:) 281 IF( MOD(mig( 1), 3) == 2 ) p3d( 2,:,:) = p3d( 1,:,:) 282 ENDIF 283 IF( nbondi /= 1 ) THEN 284 IF( MOD(mig(jpi-2), 3) == 1 ) p3d( jpi,:,:) = p3d(jpi-1,:,:) 285 IF( MOD(mig(jpi-2), 3) == 0 ) p3d(jpi-1,:,:) = p3d( jpi,:,:) 286 ENDIF 287 IF( nbondj /= -1 ) THEN 288 IF( MOD(mjg( 1), 3) == 1 ) p3d(:, 1,:) = p3d(:, 2,:) 289 IF( MOD(mjg( 1), 3) == 2 ) p3d(:, 2,:) = p3d(:, 1,:) 290 ENDIF 291 IF( nbondj /= 1 ) THEN 292 IF( MOD(mjg(jpj-2), 3) == 1 ) p3d(:, jpj,:) = p3d(:,jpj-1,:) 293 IF( MOD(mjg(jpj-2), 3) == 0 ) p3d(:,jpj-1,:) = p3d(:, jpj,:) 294 ENDIF 295 CALL lbc_lnk( 'lib_fortran', p3d, 'T', 1. ) 296 297 END SUBROUTINE sum3x3_3d 298 196 299 197 300 SUBROUTINE DDPDD( ydda, yddb ) -
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/TRP/trcrad.F90
r10365 r10408 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( 'trcrad', MIN( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:) ) 153 ztrcorn = glob_sum( 'trcrad', MIN( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:) ) 154 ! 155 IF( ztrcorb /= 0 ) THEN 156 ztrmasb = glob_sum( 'trcrad', MAX( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:) ) 157 zcoef = 1. + ztrcorb / ztrmasb 158 DO jk = 1, jpkm1 159 ptrb(:,:,jk,jn) = MAX( 0., ptrb(:,:,jk,jn) ) 160 ptrb(:,:,jk,jn) = ptrb(:,:,jk,jn) * zcoef * tmask(:,:,jk) 161 END DO 162 ENDIF 163 ! 164 IF( ztrcorn /= 0 ) THEN 165 ztrmasn = glob_sum( 'trcrad', MAX( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:) ) 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 CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb ) ! Asselin-like trend handling 179 CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn ) ! standard trend handling 180 ! 181 ENDIF 182 ! 183 END DO 184 ! 185 ELSE !== total CFC content is NOT strictly preserved ==! 186 ! 187 DO jn = jp_sms0, jp_sms1 188 ! 189 IF( l_trdtrc ) THEN 190 ztrtrdb(:,:,:) = ptrb(:,:,:,jn) ! save input trb for trend computation 191 ztrtrdn(:,:,:) = ptrn(:,:,:,jn) ! save input trn for trend computation 192 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 193 161 ! 194 162 DO jk = 1, jpkm1 195 163 DO jj = 1, jpj 196 164 DO ji = 1, jpi 197 ptrn(ji,jj,jk,jn) = MAX( 0. , ptrn(ji,jj,jk,jn) ) 198 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 199 179 END DO 200 180 END DO … … 202 182 ! 203 183 IF( l_trdtrc ) THEN 204 ! 205 zs2rdt = 1. / ( 2. * rdt * REAL( nn_dttrc, wp ) ) 206 ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt 207 ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt 208 CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb ) ! Asselin-like trend handling 209 CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn ) ! standard trend handling 210 ! 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 211 273 ENDIF 212 274 ! … … 215 277 ENDIF 216 278 ! 217 IF( l_trdtrc ) DEALLOCATE( ztrtrd b, ztrtrdn)279 IF( l_trdtrc ) DEALLOCATE( ztrtrd ) 218 280 ! 219 281 END SUBROUTINE trc_rad_sms
Note: See TracChangeset
for help on using the changeset viewer.
