- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/TOP/TRP/trcsink.F90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/TOP/TRP/trcsink.F90
r12178 r12928 24 24 INTEGER, PUBLIC :: nitermax !: Maximum number of iterations for sinking 25 25 26 !! * Substitutions 27 # include "do_loop_substitute.h90" 26 28 !!---------------------------------------------------------------------- 27 29 !! NEMO/TOP 4.0 , NEMO Consortium (2018) … … 35 37 !!---------------------------------------------------------------------- 36 38 37 SUBROUTINE trc_sink ( kt, pwsink, psinkflx, jp_tra, rsfact )39 SUBROUTINE trc_sink ( kt, Kbb, Kmm, pwsink, psinkflx, jp_tra, rsfact ) 38 40 !!--------------------------------------------------------------------- 39 41 !! *** ROUTINE trc_sink *** … … 45 47 !!--------------------------------------------------------------------- 46 48 INTEGER , INTENT(in) :: kt 49 INTEGER , INTENT(in) :: Kbb, Kmm 47 50 INTEGER , INTENT(in) :: jp_tra ! tracer index index 48 51 REAL(wp), INTENT(in) :: rsfact ! time step duration … … 70 73 iiter(:,:) = 1 71 74 ELSE 72 DO jj = 1, jpj 73 DO ji = 1, jpi 74 iiter(ji,jj) = 1 75 DO jk = 1, jpkm1 76 IF( tmask(ji,jj,jk) == 1.0 ) THEN 77 zwsmax = 0.5 * e3t_n(ji,jj,jk) * rday / rsfact 78 iiter(ji,jj) = MAX( iiter(ji,jj), INT( pwsink(ji,jj,jk) / zwsmax ) ) 79 ENDIF 80 END DO 81 END DO 82 END DO 75 DO_2D_11_11 76 iiter(ji,jj) = 1 77 DO jk = 1, jpkm1 78 IF( tmask(ji,jj,jk) == 1.0 ) THEN 79 zwsmax = 0.5 * e3t(ji,jj,jk,Kmm) * rday / rsfact 80 iiter(ji,jj) = MAX( iiter(ji,jj), INT( pwsink(ji,jj,jk) / zwsmax ) ) 81 ENDIF 82 END DO 83 END_2D 83 84 iiter(:,:) = MIN( iiter(:,:), nitermax ) 84 85 ENDIF 85 86 86 DO jk = 1,jpkm1 87 DO jj = 1, jpj 88 DO ji = 1, jpi 89 IF( tmask(ji,jj,jk) == 1.0 ) THEN 90 zwsmax = 0.5 * e3t_n(ji,jj,jk) * rday / rsfact 91 zwsink(ji,jj,jk) = MIN( pwsink(ji,jj,jk), zwsmax * REAL( iiter(ji,jj), wp ) ) 92 ELSE 93 ! provide a default value so there is no use of undefinite value in trc_sink2 for zwsink2 initialization 94 zwsink(ji,jj,jk) = 0. 95 ENDIF 96 END DO 97 END DO 98 END DO 87 DO_3D_11_11( 1,jpkm1 ) 88 IF( tmask(ji,jj,jk) == 1.0 ) THEN 89 zwsmax = 0.5 * e3t(ji,jj,jk,Kmm) * rday / rsfact 90 zwsink(ji,jj,jk) = MIN( pwsink(ji,jj,jk), zwsmax * REAL( iiter(ji,jj), wp ) ) 91 ELSE 92 ! provide a default value so there is no use of undefinite value in trc_sink2 for zwsink2 initialization 93 zwsink(ji,jj,jk) = 0. 94 ENDIF 95 END_3D 99 96 100 97 ! Initializa to zero all the sinking arrays … … 104 101 ! Compute the sedimentation term using trc_sink2 for the considered sinking particle 105 102 ! ----------------------------------------------------- 106 CALL trc_sink2( zwsink, psinkflx, jp_tra, iiter, rsfact )103 CALL trc_sink2( Kbb, Kmm, zwsink, psinkflx, jp_tra, iiter, rsfact ) 107 104 ! 108 105 IF( ln_timing ) CALL timing_stop('trc_sink') … … 110 107 END SUBROUTINE trc_sink 111 108 112 SUBROUTINE trc_sink2( pwsink, psinkflx, jp_tra, kiter, rsfact )109 SUBROUTINE trc_sink2( Kbb, Kmm, pwsink, psinkflx, jp_tra, kiter, rsfact ) 113 110 !!--------------------------------------------------------------------- 114 111 !! *** ROUTINE trc_sink2 *** … … 121 118 !! transport term, i.e. div(u*tra). 122 119 !!--------------------------------------------------------------------- 120 INTEGER, INTENT(in ) :: Kbb, Kmm ! time level indices 123 121 INTEGER, INTENT(in ) :: jp_tra ! tracer index index 124 122 REAL(wp), INTENT(in ) :: rsfact ! duration of time step … … 136 134 ztraz(:,:,:) = 0.e0 137 135 zakz (:,:,:) = 0.e0 138 ztrb (:,:,:) = tr b(:,:,:,jp_tra)136 ztrb (:,:,:) = tr(:,:,:,jp_tra,Kbb) 139 137 140 138 DO jk = 1, jpkm1 … … 147 145 DO jn = 1, 2 148 146 ! first guess of the slopes interior values 149 DO jj = 1, jpj 150 DO ji = 1, jpi 151 ! 152 zstep = rsfact / REAL( kiter(ji,jj), wp ) / 2. 153 ! 154 DO jk = 2, jpkm1 155 ztraz(ji,jj,jk) = ( trb(ji,jj,jk-1,jp_tra) - trb(ji,jj,jk,jp_tra) ) * tmask(ji,jj,jk) 156 END DO 157 ztraz(ji,jj,1 ) = 0.0 158 ztraz(ji,jj,jpk) = 0.0 159 160 ! slopes 161 DO jk = 2, jpkm1 162 zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) ) 163 zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign 164 END DO 165 166 ! Slopes limitation 167 DO jk = 2, jpkm1 168 zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) * & 169 & MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) ) 170 END DO 171 172 ! vertical advective flux 173 DO jk = 1, jpkm1 174 zigma = zwsink2(ji,jj,jk+1) * zstep / e3w_n(ji,jj,jk+1) 175 zew = zwsink2(ji,jj,jk+1) 176 psinkflx(ji,jj,jk+1) = -zew * ( trb(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep 177 END DO 178 ! 179 ! Boundary conditions 180 psinkflx(ji,jj,1 ) = 0.e0 181 psinkflx(ji,jj,jpk) = 0.e0 182 183 DO jk=1,jpkm1 184 zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk) 185 trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + zflx 186 END DO 187 END DO 188 END DO 147 DO_2D_11_11 148 ! 149 zstep = rsfact / REAL( kiter(ji,jj), wp ) / 2. 150 ! 151 DO jk = 2, jpkm1 152 ztraz(ji,jj,jk) = ( tr(ji,jj,jk-1,jp_tra,Kbb) - tr(ji,jj,jk,jp_tra,Kbb) ) * tmask(ji,jj,jk) 153 END DO 154 ztraz(ji,jj,1 ) = 0.0 155 ztraz(ji,jj,jpk) = 0.0 156 157 ! slopes 158 DO jk = 2, jpkm1 159 zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) ) 160 zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign 161 END DO 162 163 ! Slopes limitation 164 DO jk = 2, jpkm1 165 zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) * & 166 & MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) ) 167 END DO 168 169 ! vertical advective flux 170 DO jk = 1, jpkm1 171 zigma = zwsink2(ji,jj,jk+1) * zstep / e3w(ji,jj,jk+1,Kmm) 172 zew = zwsink2(ji,jj,jk+1) 173 psinkflx(ji,jj,jk+1) = -zew * ( tr(ji,jj,jk,jp_tra,Kbb) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep 174 END DO 175 ! 176 ! Boundary conditions 177 psinkflx(ji,jj,1 ) = 0.e0 178 psinkflx(ji,jj,jpk) = 0.e0 179 180 DO jk=1,jpkm1 181 zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm) 182 tr(ji,jj,jk,jp_tra,Kbb) = tr(ji,jj,jk,jp_tra,Kbb) + zflx 183 END DO 184 END_2D 189 185 END DO 190 186 191 DO jk = 1,jpkm1 192 DO jj = 1,jpj 193 DO ji = 1, jpi 194 zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk) 195 ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx 196 END DO 197 END DO 198 END DO 199 200 trb(:,:,:,jp_tra) = ztrb(:,:,:) 187 DO_3D_11_11( 1,jpkm1 ) 188 zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm) 189 ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx 190 END_3D 191 192 tr(:,:,:,jp_tra,Kbb) = ztrb(:,:,:) 201 193 psinkflx(:,:,:) = 2. * psinkflx(:,:,:) 202 194 ! … … 216 208 !!---------------------------------------------------------------------- 217 209 ! 218 REWIND( numnat_ref ) ! namtrc_rad in reference namelist219 210 READ ( numnat_ref, namtrc_snk, IOSTAT = ios, ERR = 907) 220 211 907 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_snk in reference namelist' ) 221 REWIND( numnat_cfg ) ! namtrc_rad in configuration namelist222 212 READ ( numnat_cfg, namtrc_snk, IOSTAT = ios, ERR = 908 ) 223 213 908 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtrc_snk in configuration namelist' )
Note: See TracChangeset
for help on using the changeset viewer.
