Changeset 10375 for NEMO/trunk/src/TOP/PISCES/P4Z/p4zsink.F90
- Timestamp:
- 2018-12-07T17:27:26+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/TOP/PISCES/P4Z/p4zsink.F90
r10362 r10375 16 16 USE trc ! passive tracers common variables 17 17 USE sms_pisces ! PISCES Source Minus Sink variables 18 USE trcsink ! General routine to compute sedimentation 18 19 USE prtctl_trc ! print control for debugging 19 20 USE iom ! I/O manager … … 59 60 !!--------------------------------------------------------------------- 60 61 INTEGER, INTENT(in) :: kt, knt 61 INTEGER :: ji, jj, jk, jit 62 INTEGER :: iiter1, iiter2 63 REAL(wp) :: zagg1, zagg2, zagg3, zagg4 64 REAL(wp) :: zagg , zaggfe, zaggdoc, zaggdoc2, zaggdoc3 65 REAL(wp) :: zfact, zwsmax, zmax 62 INTEGER :: ji, jj, jk 66 63 CHARACTER (len=25) :: charout 64 REAL(wp) :: zmax, zfact 67 65 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d 68 66 REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: zw2d … … 70 68 ! 71 69 IF( ln_timing ) CALL timing_start('p4z_sink') 72 73 70 74 71 ! Initialization of some global variables … … 97 94 98 95 ! 99 ! OA This is (I hope) a temporary solution for the problem that may100 ! OA arise in specific situation where the CFL criterion is broken101 ! OA for vertical sedimentation of particles. To avoid this, a time102 ! OA splitting algorithm has been coded. A specific maximum103 ! OA iteration number is provided and may be specified in the namelist104 ! OA This is to avoid very large iteration number when explicit free105 ! OA surface is used (for instance). When niter?max is set to 1,106 ! OA this computation is skipped. The crude old threshold method is107 ! OA then applied. This also happens when niter exceeds nitermax.108 IF( MAX( niter1max, niter2max ) == 1 ) THEN109 iiter1 = 1110 iiter2 = 1111 ELSE112 iiter1 = 1113 iiter2 = 1114 DO jk = 1, jpkm1115 DO jj = 1, jpj116 DO ji = 1, jpi117 IF( tmask(ji,jj,jk) == 1) THEN118 zwsmax = 0.5 * e3t_n(ji,jj,jk) / xstep119 iiter1 = MAX( iiter1, INT( wsbio3(ji,jj,jk) / zwsmax ) )120 iiter2 = MAX( iiter2, INT( wsbio4(ji,jj,jk) / zwsmax ) )121 ENDIF122 END DO123 END DO124 END DO125 IF( lk_mpp ) THEN126 CALL mpp_max( iiter1 )127 CALL mpp_max( iiter2 )128 ENDIF129 iiter1 = MIN( iiter1, niter1max )130 iiter2 = MIN( iiter2, niter2max )131 ENDIF132 133 DO jk = 1,jpkm1134 DO jj = 1, jpj135 DO ji = 1, jpi136 IF( tmask(ji,jj,jk) == 1 ) THEN137 zwsmax = 0.5 * e3t_n(ji,jj,jk) / xstep138 wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax * REAL( iiter1, wp ) )139 wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax * REAL( iiter2, wp ) )140 ENDIF141 END DO142 END DO143 END DO144 145 96 ! Initializa to zero all the sinking arrays 146 97 ! ----------------------------------------- … … 154 105 ! Compute the sedimentation term using p4zsink2 for all the sinking particles 155 106 ! ----------------------------------------------------- 156 DO jit = 1, iiter1 157 CALL p4z_sink2( wsbio3, sinking , jppoc, iiter1 ) 158 CALL p4z_sink2( wsbio3, sinkfer , jpsfe, iiter1 ) 159 END DO 160 161 DO jit = 1, iiter2 162 CALL p4z_sink2( wsbio4, sinking2, jpgoc, iiter2 ) 163 CALL p4z_sink2( wsbio4, sinkfer2, jpbfe, iiter2 ) 164 CALL p4z_sink2( wsbio4, sinksil , jpgsi, iiter2 ) 165 CALL p4z_sink2( wsbio4, sinkcal , jpcal, iiter2 ) 166 END DO 107 CALL trc_sink( kt, wsbio3, sinking , jppoc, rfact2 ) 108 CALL trc_sink( kt, wsbio3, sinkfer , jpsfe, rfact2 ) 109 CALL trc_sink( kt, wsbio4, sinking2, jpgoc, rfact2 ) 110 CALL trc_sink( kt, wsbio4, sinkfer2, jpbfe, rfact2 ) 111 CALL trc_sink( kt, wsbio4, sinksil , jpgsi, rfact2 ) 112 CALL trc_sink( kt, wsbio4, sinkcal , jpcal, rfact2 ) 167 113 168 114 IF( ln_p5z ) THEN … … 174 120 ! Compute the sedimentation term using p4zsink2 for all the sinking particles 175 121 ! ----------------------------------------------------- 176 DO jit = 1, iiter1 177 CALL p4z_sink2( wsbio3, sinkingn , jppon, iiter1 ) 178 CALL p4z_sink2( wsbio3, sinkingp , jppop, iiter1 ) 179 END DO 180 181 DO jit = 1, iiter2 182 CALL p4z_sink2( wsbio4, sinking2n, jpgon, iiter2 ) 183 CALL p4z_sink2( wsbio4, sinking2p, jpgop, iiter2 ) 184 END DO 122 CALL trc_sink( kt, wsbio3, sinkingn , jppon, rfact2 ) 123 CALL trc_sink( kt, wsbio3, sinkingp , jppop, rfact2 ) 124 CALL trc_sink( kt, wsbio4, sinking2n, jpgon, rfact2 ) 125 CALL trc_sink( kt, wsbio4, sinking2p, jpgop, rfact2 ) 185 126 ENDIF 186 127 187 128 IF( ln_ligand ) THEN 188 129 wsfep (:,:,:) = wfep 189 DO jk = 1,jpkm1190 DO jj = 1, jpj191 DO ji = 1, jpi192 IF( tmask(ji,jj,jk) == 1 ) THEN193 zwsmax = 0.5 * e3t_n(ji,jj,jk) / xstep194 wsfep(ji,jj,jk) = MIN( wsfep(ji,jj,jk), zwsmax * REAL( iiter1, wp ) )195 ENDIF196 END DO197 END DO198 END DO199 130 ! 200 131 sinkfep(:,:,:) = 0.e0 201 DO jit = 1, iiter1 202 CALL p4z_sink2( wsfep, sinkfep , jpfep, iiter1 ) 203 END DO 132 CALL trc_sink( kt, wsfep, sinkfep , jpfep, rfact2 ) 204 133 ENDIF 205 134 … … 281 210 END SUBROUTINE p4z_sink_init 282 211 283 284 SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra, kiter )285 !!---------------------------------------------------------------------286 !! *** ROUTINE p4z_sink2 ***287 !!288 !! ** Purpose : Compute the sedimentation terms for the various sinking289 !! particles. The scheme used to compute the trends is based290 !! on MUSCL.291 !!292 !! ** Method : - this ROUTINE compute not exactly the advection but the293 !! transport term, i.e. div(u*tra).294 !!---------------------------------------------------------------------295 INTEGER , INTENT(in ) :: jp_tra ! tracer index index296 INTEGER , INTENT(in ) :: kiter ! number of iterations for time-splitting297 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwsink ! sinking speed298 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: psinkflx ! sinking fluxe299 !300 INTEGER :: ji, jj, jk, jn301 REAL(wp) :: zigma,zew,zign, zflx, zstep302 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztraz, zakz, zwsink2, ztrb303 !!---------------------------------------------------------------------304 !305 IF( ln_timing ) CALL timing_start('p4z_sink2')306 !307 zstep = rfact2 / REAL( kiter, wp ) / 2.308 309 ztraz(:,:,:) = 0.e0310 zakz (:,:,:) = 0.e0311 ztrb (:,:,:) = trb(:,:,:,jp_tra)312 313 DO jk = 1, jpkm1314 zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1)315 END DO316 zwsink2(:,:,1) = 0.e0317 318 319 ! Vertical advective flux320 DO jn = 1, 2321 ! first guess of the slopes interior values322 DO jk = 2, jpkm1323 ztraz(:,:,jk) = ( trb(:,:,jk-1,jp_tra) - trb(:,:,jk,jp_tra) ) * tmask(:,:,jk)324 END DO325 ztraz(:,:,1 ) = 0.0326 ztraz(:,:,jpk) = 0.0327 328 ! slopes329 DO jk = 2, jpkm1330 DO jj = 1,jpj331 DO ji = 1, jpi332 zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )333 zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign334 END DO335 END DO336 END DO337 338 ! Slopes limitation339 DO jk = 2, jpkm1340 DO jj = 1, jpj341 DO ji = 1, jpi342 zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) * &343 & MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )344 END DO345 END DO346 END DO347 348 ! vertical advective flux349 DO jk = 1, jpkm1350 DO jj = 1, jpj351 DO ji = 1, jpi352 zigma = zwsink2(ji,jj,jk+1) * zstep / e3w_n(ji,jj,jk+1)353 zew = zwsink2(ji,jj,jk+1)354 psinkflx(ji,jj,jk+1) = -zew * ( trb(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep355 END DO356 END DO357 END DO358 !359 ! Boundary conditions360 psinkflx(:,:,1 ) = 0.e0361 psinkflx(:,:,jpk) = 0.e0362 363 DO jk=1,jpkm1364 DO jj = 1,jpj365 DO ji = 1, jpi366 zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)367 trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + zflx368 END DO369 END DO370 END DO371 372 ENDDO373 374 DO jk = 1,jpkm1375 DO jj = 1,jpj376 DO ji = 1, jpi377 zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)378 ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx379 END DO380 END DO381 END DO382 383 trb(:,:,:,jp_tra) = ztrb(:,:,:)384 psinkflx(:,:,:) = 2. * psinkflx(:,:,:)385 !386 IF( ln_timing ) CALL timing_stop('p4z_sink2')387 !388 END SUBROUTINE p4z_sink2389 390 391 212 INTEGER FUNCTION p4z_sink_alloc() 392 213 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.