![(please configure the [header_logo] section in trac.ini)](/dynamico/chrome/site/your_project_logo.png){kind=logo}
Changeset 174 for codes/icosagcm/trunk/src/advect.f90
- Timestamp:
- 10/07/13 17:48:24 (11 years ago)
- File:
-
- 1 edited
-
codes/icosagcm/trunk/src/advect.f90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/advect.f90
r151 r174 17 17 REAL(rstd),INTENT(OUT) :: one_over_sqrt_leng(iim*jjm) 18 18 19 INTEGER :: i,j,n 20 21 DO j=jj_begin-1,jj_end+1 22 DO i=ii_begin-1,ii_end+1 23 n=(j-1)*iim+i 24 25 CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),normal(n+u_right,:)) 26 normal(n+u_right,:)=normal(n+u_right,:)/sqrt(sum(normal(n+u_right,:)**2)+1e-50)*ne(n,right) 27 28 CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),normal(n+u_lup,:)) 29 normal(n+u_lup,:)=normal(n+u_lup,:)/sqrt(sum(normal(n+u_lup,:)**2)+1e-50)*ne(n,lup) 30 31 CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),normal(n+u_ldown,:)) 32 normal(n+u_ldown,:)=normal(n+u_ldown,:)/sqrt(sum(normal(n+u_ldown,:)**2)+1e-50)*ne(n,ldown) 33 34 tangent(n+u_right,:)=xyz_v(n+z_rup,:)-xyz_v(n+z_rdown,:) 35 tangent(n+u_right,:)=tangent(n+u_right,:)/sqrt(sum(tangent(n+u_right,:)**2)+1e-50) 36 37 tangent(n+u_lup,:)=xyz_v(n+z_lup,:)-xyz_v(n+z_up,:) 38 tangent(n+u_lup,:)=tangent(n+u_lup,:)/sqrt(sum(tangent(n+u_lup,:)**2)+1e-50) 39 40 tangent(n+u_ldown,:)=xyz_v(n+z_down,:)-xyz_v(n+z_ldown,:) 41 tangent(n+u_ldown,:)=tangent(n+u_ldown,:)/sqrt(sum(tangent(n+u_ldown,:)**2)+1e-50) 42 43 one_over_sqrt_leng(n) = 1./sqrt(max(sum((xyz_v(n+z_up,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_down,:) - xyz_i(n,:))**2), & 44 sum((xyz_v(n+z_rup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_rdown,:) - xyz_i(n,:))**2), & 45 sum((xyz_v(n+z_lup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_ldown,:) - xyz_i(n,:))**2)) ) 46 ENDDO 19 INTEGER :: ij 20 21 !$SIMD 22 DO ij=ij_begin,ij_end 23 24 CALL cross_product2(xyz_v(ij+z_rdown,:),xyz_v(ij+z_rup,:),normal(ij+u_right,:)) 25 normal(ij+u_right,:)=normal(ij+u_right,:)/sqrt(sum(normal(ij+u_right,:)**2)+1e-50)*ne(ij,right) 26 27 CALL cross_product2(xyz_v(ij+z_up,:),xyz_v(ij+z_lup,:),normal(ij+u_lup,:)) 28 normal(ij+u_lup,:)=normal(ij+u_lup,:)/sqrt(sum(normal(ij+u_lup,:)**2)+1e-50)*ne(ij,lup) 29 30 CALL cross_product2(xyz_v(ij+z_ldown,:),xyz_v(ij+z_down,:),normal(ij+u_ldown,:)) 31 normal(ij+u_ldown,:)=normal(ij+u_ldown,:)/sqrt(sum(normal(ij+u_ldown,:)**2)+1e-50)*ne(ij,ldown) 32 33 tangent(ij+u_right,:)=xyz_v(ij+z_rup,:)-xyz_v(ij+z_rdown,:) 34 tangent(ij+u_right,:)=tangent(ij+u_right,:)/sqrt(sum(tangent(ij+u_right,:)**2)+1e-50) 35 36 tangent(ij+u_lup,:)=xyz_v(ij+z_lup,:)-xyz_v(ij+z_up,:) 37 tangent(ij+u_lup,:)=tangent(ij+u_lup,:)/sqrt(sum(tangent(ij+u_lup,:)**2)+1e-50) 38 39 tangent(ij+u_ldown,:)=xyz_v(ij+z_down,:)-xyz_v(ij+z_ldown,:) 40 tangent(ij+u_ldown,:)=tangent(ij+u_ldown,:)/sqrt(sum(tangent(ij+u_ldown,:)**2)+1e-50) 41 42 one_over_sqrt_leng(ij) = 1./sqrt(max(sum((xyz_v(ij+z_up,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_down,:) - xyz_i(ij,:))**2), & 43 sum((xyz_v(ij+z_rup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_rdown,:) - xyz_i(ij,:))**2), & 44 sum((xyz_v(ij+z_lup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_ldown,:) - xyz_i(ij,:))**2)) ) 47 45 ENDDO 48 46 … … 64 62 REAL(rstd) :: ar(iim*jjm) 65 63 REAL(rstd) :: gradtri(2*iim*jjm,llm,3) 66 INTEGER :: i ,j,n,k,ind,l64 INTEGER :: ij,k,ind,l 67 65 68 66 CALL trace_start("compute_gradq3d") … … 74 72 ! arr = area of triangle joining centroids of hexagons 75 73 DO l = ll_begin,ll_end 76 DO j=jj_begin-1,jj_end+1 77 DO i=ii_begin-1,ii_end+1 78 n=(j-1)*iim+i 79 CALL gradq(n,l,n+t_rup,n+t_lup,n+z_up,qi,gradtri(n+z_up,l,:),arr(n+z_up)) 80 CALL gradq(n,l,n+t_ldown,n+t_rdown,n+z_down,qi,gradtri(n+z_down,l,:),arr(n+z_down)) 81 END DO 74 !$SIMD 75 DO ij=ij_begin_ext,ij_end_ext 76 CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,:),arr(ij+z_up)) 77 CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,:),arr(ij+z_down)) 82 78 END DO 83 79 END DO 84 80 85 DO j=jj_begin-1,jj_end+1 86 DO i=ii_begin-1,ii_end+1 87 n=(j-1)*iim+i 88 ar(n) = arr(n+z_up)+arr(n+z_lup)+arr(n+z_ldown)+arr(n+z_down)+arr(n+z_rdown)+arr(n+z_rup)+1.e-50 89 ENDDO 81 !$SIMD 82 DO ij=ij_begin,ij_end 83 ar(ij) = arr(ij+z_up)+arr(ij+z_lup)+arr(ij+z_ldown)+arr(ij+z_down)+arr(ij+z_rdown)+arr(ij+z_rup)+1.e-50 90 84 ENDDO 91 85 92 86 DO k=1,3 93 87 DO l =ll_begin,ll_end 94 DO j=jj_begin,jj_end 95 DO i=ii_begin,ii_end 96 n=(j-1)*iim+i 97 gradq3d(n,l,k) = ( gradtri(n+z_up,l,k) + gradtri(n+z_down,l,k) + & 98 gradtri(n+z_rup,l,k) + gradtri(n+z_ldown,l,k) + & 99 gradtri(n+z_lup,l,k)+ gradtri(n+z_rdown,l,k) ) / ar(n) 100 END DO 88 !$SIMD 89 DO ij=ij_begin,ij_end 90 gradq3d(ij,l,k) = ( gradtri(ij+z_up,l,k) + gradtri(ij+z_down,l,k) + & 91 gradtri(ij+z_rup,l,k) + gradtri(ij+z_ldown,l,k) + & 92 gradtri(ij+z_lup,l,k)+ gradtri(ij+z_rdown,l,k) ) / ar(ij) 101 93 END DO 102 94 END DO … … 105 97 !============================================================================================= LIMITING 106 98 DO l =ll_begin,ll_end 107 DO j=jj_begin,jj_end 108 DO i=ii_begin,ii_end 109 n=(j-1)*iim+i 110 maggrd = dot_product(gradq3d(n,l,:),gradq3d(n,l,:)) 99 !$SIMD 100 DO ij=ij_begin,ij_end 101 maggrd = dot_product(gradq3d(ij,l,:),gradq3d(ij,l,:)) 111 102 maggrd = sqrt(maggrd) 112 maxq_c = qi( n,l) + maggrd*one_over_sqrt_leng(n)113 minq_c = qi( n,l) - maggrd*one_over_sqrt_leng(n)114 maxq = max(qi( n,l),qi(n+t_right,l),qi(n+t_lup,l),qi(n+t_rup,l),qi(n+t_left,l), &115 qi( n+t_rdown,l),qi(n+t_ldown,l))116 minq = min(qi( n,l),qi(n+t_right,l),qi(n+t_lup,l),qi(n+t_rup,l),qi(n+t_left,l), &117 qi( n+t_rdown,l),qi(n+t_ldown,l))118 alphamx = (maxq - qi( n,l)) ; alphamx = alphamx/(maxq_c - qi(n,l) )103 maxq_c = qi(ij,l) + maggrd*one_over_sqrt_leng(ij) 104 minq_c = qi(ij,l) - maggrd*one_over_sqrt_leng(ij) 105 maxq = max(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), & 106 qi(ij+t_rdown,l),qi(ij+t_ldown,l)) 107 minq = min(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), & 108 qi(ij+t_rdown,l),qi(ij+t_ldown,l)) 109 alphamx = (maxq - qi(ij,l)) ; alphamx = alphamx/(maxq_c - qi(ij,l) ) 119 110 alphamx = max(alphamx,0.0) 120 alphami = (minq - qi( n,l)); alphami = alphami/(minq_c - qi(n,l))111 alphami = (minq - qi(ij,l)); alphami = alphami/(minq_c - qi(ij,l)) 121 112 alphami = max(alphami,0.0) 122 113 alpha = min(alphamx,alphami,1.0) 123 gradq3d(n,l,:) = alpha*gradq3d(n,l,:) 124 END DO 114 gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:) 125 115 END DO 126 116 END DO … … 141 131 142 132 REAL(rstd) :: v_e(3), up_e, qe, ed(3) 143 INTEGER :: i ,j,n,l133 INTEGER :: ij,l 144 134 145 135 CALL trace_start("compute_backward_traj") … … 149 139 ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau 150 140 DO l = ll_begin,ll_end 151 DO j=jj_begin-1,jj_end+1 152 DO i=ii_begin-1,ii_end+1 153 n=(j-1)*iim+i 154 155 up_e =1/de(n+u_right)*( & 156 wee(n+u_right,1,1)*le(n+u_rup)*ue(n+u_rup,l)+ & 157 wee(n+u_right,2,1)*le(n+u_lup)*ue(n+u_lup,l)+ & 158 wee(n+u_right,3,1)*le(n+u_left)*ue(n+u_left,l)+ & 159 wee(n+u_right,4,1)*le(n+u_ldown)*ue(n+u_ldown,l)+ & 160 wee(n+u_right,5,1)*le(n+u_rdown)*ue(n+u_rdown,l)+ & 161 wee(n+u_right,1,2)*le(n+t_right+u_ldown)*ue(n+t_right+u_ldown,l)+ & 162 wee(n+u_right,2,2)*le(n+t_right+u_rdown)*ue(n+t_right+u_rdown,l)+ & 163 wee(n+u_right,3,2)*le(n+t_right+u_right)*ue(n+t_right+u_right,l)+ & 164 wee(n+u_right,4,2)*le(n+t_right+u_rup)*ue(n+t_right+u_rup,l)+ & 165 wee(n+u_right,5,2)*le(n+t_right+u_lup)*ue(n+t_right+u_lup,l) & 141 !$SIMD 142 DO ij=ij_begin,ij_end 143 up_e =1/de(ij+u_right)*( & 144 wee(ij+u_right,1,1)*le(ij+u_rup)*ue(ij+u_rup,l)+ & 145 wee(ij+u_right,2,1)*le(ij+u_lup)*ue(ij+u_lup,l)+ & 146 wee(ij+u_right,3,1)*le(ij+u_left)*ue(ij+u_left,l)+ & 147 wee(ij+u_right,4,1)*le(ij+u_ldown)*ue(ij+u_ldown,l)+ & 148 wee(ij+u_right,5,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+ & 149 wee(ij+u_right,1,2)*le(ij+t_right+u_ldown)*ue(ij+t_right+u_ldown,l)+ & 150 wee(ij+u_right,2,2)*le(ij+t_right+u_rdown)*ue(ij+t_right+u_rdown,l)+ & 151 wee(ij+u_right,3,2)*le(ij+t_right+u_right)*ue(ij+t_right+u_right,l)+ & 152 wee(ij+u_right,4,2)*le(ij+t_right+u_rup)*ue(ij+t_right+u_rup,l)+ & 153 wee(ij+u_right,5,2)*le(ij+t_right+u_lup)*ue(ij+t_right+u_lup,l) & 166 154 ) 167 v_e = ue( n+u_right,l)*normal(n+u_right,:) + up_e*tangent(n+u_right,:)168 cc( n+u_right,l,:) = xyz_e(n+u_right,:) - v_e*tau169 170 up_e=1/de( n+u_lup)*( &171 wee( n+u_lup,1,1)*le(n+u_left)*ue(n+u_left,l)+ &172 wee( n+u_lup,2,1)*le(n+u_ldown)*ue(n+u_ldown,l)+ &173 wee( n+u_lup,3,1)*le(n+u_rdown)*ue(n+u_rdown,l)+ &174 wee( n+u_lup,4,1)*le(n+u_right)*ue(n+u_right,l)+ &175 wee( n+u_lup,5,1)*le(n+u_rup)*ue(n+u_rup,l)+ &176 wee( n+u_lup,1,2)*le(n+t_lup+u_right)*ue(n+t_lup+u_right,l)+ &177 wee( n+u_lup,2,2)*le(n+t_lup+u_rup)*ue(n+t_lup+u_rup,l)+ &178 wee( n+u_lup,3,2)*le(n+t_lup+u_lup)*ue(n+t_lup+u_lup,l)+ &179 wee( n+u_lup,4,2)*le(n+t_lup+u_left)*ue(n+t_lup+u_left,l)+ &180 wee( n+u_lup,5,2)*le(n+t_lup+u_ldown)*ue(n+t_lup+u_ldown,l) &155 v_e = ue(ij+u_right,l)*normal(ij+u_right,:) + up_e*tangent(ij+u_right,:) 156 cc(ij+u_right,l,:) = xyz_e(ij+u_right,:) - v_e*tau 157 158 up_e=1/de(ij+u_lup)*( & 159 wee(ij+u_lup,1,1)*le(ij+u_left)*ue(ij+u_left,l)+ & 160 wee(ij+u_lup,2,1)*le(ij+u_ldown)*ue(ij+u_ldown,l)+ & 161 wee(ij+u_lup,3,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+ & 162 wee(ij+u_lup,4,1)*le(ij+u_right)*ue(ij+u_right,l)+ & 163 wee(ij+u_lup,5,1)*le(ij+u_rup)*ue(ij+u_rup,l)+ & 164 wee(ij+u_lup,1,2)*le(ij+t_lup+u_right)*ue(ij+t_lup+u_right,l)+ & 165 wee(ij+u_lup,2,2)*le(ij+t_lup+u_rup)*ue(ij+t_lup+u_rup,l)+ & 166 wee(ij+u_lup,3,2)*le(ij+t_lup+u_lup)*ue(ij+t_lup+u_lup,l)+ & 167 wee(ij+u_lup,4,2)*le(ij+t_lup+u_left)*ue(ij+t_lup+u_left,l)+ & 168 wee(ij+u_lup,5,2)*le(ij+t_lup+u_ldown)*ue(ij+t_lup+u_ldown,l) & 181 169 ) 182 v_e = ue(n+u_lup,l)*normal(n+u_lup,:) + up_e*tangent(n+u_lup,:) 183 cc(n+u_lup,l,:) = xyz_e(n+u_lup,:) - v_e*tau 184 185 up_e=1/de(n+u_ldown)*( & 186 wee(n+u_ldown,1,1)*le(n+u_rdown)*ue(n+u_rdown,l)+ & 187 wee(n+u_ldown,2,1)*le(n+u_right)*ue(n+u_right,l)+ & 188 wee(n+u_ldown,3,1)*le(n+u_rup)*ue(n+u_rup,l)+ & 189 wee(n+u_ldown,4,1)*le(n+u_lup)*ue(n+u_lup,l)+ & 190 wee(n+u_ldown,5,1)*le(n+u_left)*ue(n+u_left,l)+ & 191 wee(n+u_ldown,1,2)*le(n+t_ldown+u_lup)*ue(n+t_ldown+u_lup,l)+ & 192 wee(n+u_ldown,2,2)*le(n+t_ldown+u_left)*ue(n+t_ldown+u_left,l)+ & 193 wee(n+u_ldown,3,2)*le(n+t_ldown+u_ldown)*ue(n+t_ldown+u_ldown,l)+ & 194 wee(n+u_ldown,4,2)*le(n+t_ldown+u_rdown)*ue(n+t_ldown+u_rdown,l)+ & 195 wee(n+u_ldown,5,2)*le(n+t_ldown+u_right)*ue(n+t_ldown+u_right,l) & 170 v_e = ue(ij+u_lup,l)*normal(ij+u_lup,:) + up_e*tangent(ij+u_lup,:) 171 cc(ij+u_lup,l,:) = xyz_e(ij+u_lup,:) - v_e*tau 172 173 174 up_e=1/de(ij+u_ldown)*( & 175 wee(ij+u_ldown,1,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+ & 176 wee(ij+u_ldown,2,1)*le(ij+u_right)*ue(ij+u_right,l)+ & 177 wee(ij+u_ldown,3,1)*le(ij+u_rup)*ue(ij+u_rup,l)+ & 178 wee(ij+u_ldown,4,1)*le(ij+u_lup)*ue(ij+u_lup,l)+ & 179 wee(ij+u_ldown,5,1)*le(ij+u_left)*ue(ij+u_left,l)+ & 180 wee(ij+u_ldown,1,2)*le(ij+t_ldown+u_lup)*ue(ij+t_ldown+u_lup,l)+ & 181 wee(ij+u_ldown,2,2)*le(ij+t_ldown+u_left)*ue(ij+t_ldown+u_left,l)+ & 182 wee(ij+u_ldown,3,2)*le(ij+t_ldown+u_ldown)*ue(ij+t_ldown+u_ldown,l)+ & 183 wee(ij+u_ldown,4,2)*le(ij+t_ldown+u_rdown)*ue(ij+t_ldown+u_rdown,l)+ & 184 wee(ij+u_ldown,5,2)*le(ij+t_ldown+u_right)*ue(ij+t_ldown+u_right,l) & 196 185 ) 197 v_e = ue(n+u_ldown,l)*normal(n+u_ldown,:) + up_e*tangent(n+u_ldown,:) 198 cc(n+u_ldown,l,:) = xyz_e(n+u_ldown,:) - v_e*tau 199 ENDDO 186 v_e = ue(ij+u_ldown,l)*normal(ij+u_ldown,:) + up_e*tangent(ij+u_ldown,:) 187 cc(ij+u_ldown,l,:) = xyz_e(ij+u_ldown,:) - v_e*tau 200 188 ENDDO 201 189 END DO … … 220 208 REAL(rstd) :: dq,dmass,qe,ed(3), newmass 221 209 REAL(rstd) :: qflux(3*iim*jjm,llm) 222 INTEGER :: i ,j,n,k,l210 INTEGER :: ij,k,l 223 211 224 212 CALL trace_start("compute_advect_horiz") … … 228 216 ! ne*hfluxt>0 iff outgoing 229 217 DO l = ll_begin,ll_end 230 DO j=jj_begin-1,jj_end+1 231 DO i=ii_begin-1,ii_end+1 232 n=(j-1)*iim+i 233 if (ne(n,right)*hfluxt(n+u_right,l)>0) then 234 ed = cc(n+u_right,l,:) - xyz_i(n,:) 235 qe = qi(n,l)+sum2(gradq3d(n,l,:),ed) 236 else 237 ed = cc(n+u_right,l,:) - xyz_i(n+t_right,:) 238 qe = qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed) 239 endif 240 qflux(n+u_right,l) = hfluxt(n+u_right,l)*qe 218 219 !$SIMD 220 DO ij=ij_begin_ext,ij_end_ext 221 222 IF (ne(ij,right)*hfluxt(ij+u_right,l)>0) THEN 223 ed = cc(ij+u_right,l,:) - xyz_i(ij,:) 224 qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 225 ELSE 226 ed = cc(ij+u_right,l,:) - xyz_i(ij+t_right,:) 227 qe = qi(ij+t_right,l)+sum2(gradq3d(ij+t_right,l,:),ed) 228 ENDIF 229 qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe 241 230 242 if (ne(n,lup)*hfluxt(n+u_lup,l)>0) then 243 ed = cc(n+u_lup,l,:) - xyz_i(n,:) 244 qe = qi(n,l)+sum2(gradq3d(n,l,:),ed) 245 else 246 ed = cc(n+u_lup,l,:) - xyz_i(n+t_lup,:) 247 qe = qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed) 248 endif 249 qflux(n+u_lup,l) = hfluxt(n+u_lup,l)*qe 250 251 if (ne(n,ldown)*hfluxt(n+u_ldown,l)>0) then 252 ed = cc(n+u_ldown,l,:) - xyz_i(n,:) 253 qe = qi(n,l)+sum2(gradq3d(n,l,:),ed) 254 else 255 ed = cc(n+u_ldown,l,:) - xyz_i(n+t_ldown,:) 256 qe = qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed) 257 endif 258 qflux(n+u_ldown,l) = hfluxt(n+u_ldown,l)*qe 259 end do 260 end do 231 IF (ne(ij,lup)*hfluxt(ij+u_lup,l)>0) THEN 232 ed = cc(ij+u_lup,l,:) - xyz_i(ij,:) 233 qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 234 ELSE 235 ed = cc(ij+u_lup,l,:) - xyz_i(ij+t_lup,:) 236 qe = qi(ij+t_lup,l)+sum2(gradq3d(ij+t_lup,l,:),ed) 237 ENDIF 238 qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe 239 240 IF (ne(ij,ldown)*hfluxt(ij+u_ldown,l)>0) THEN 241 ed = cc(ij+u_ldown,l,:) - xyz_i(ij,:) 242 qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 243 ELSE 244 ed = cc(ij+u_ldown,l,:) - xyz_i(ij+t_ldown,:) 245 qe = qi(ij+t_ldown,l)+sum2(gradq3d(ij+t_ldown,l,:),ed) 246 ENDIF 247 qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe 248 END DO 261 249 END DO 262 250 263 251 ! update q and, if update_mass, update mass 264 252 DO l = ll_begin,ll_end 265 DO j=jj_begin,jj_end 266 DO i=ii_begin,ii_end 267 n=(j-1)*iim+i 253 !$SIMD 254 DO ij=ij_begin,ij_end 268 255 ! sign convention as Ringler et al. (2010) eq. 21 269 dmass = hfluxt(n+u_right,l)*ne(n,right) & 270 + hfluxt(n+u_lup,l) *ne(n,lup) & 271 + hfluxt(n+u_ldown,l)*ne(n,ldown) & 272 + hfluxt(n+u_rup,l) *ne(n,rup) & 273 + hfluxt(n+u_left,l) *ne(n,left) & 274 + hfluxt(n+u_rdown,l)*ne(n,rdown) 275 276 dq = qflux(n+u_right,l) *ne(n,right) & 277 + qflux(n+u_lup,l) *ne(n,lup) & 278 + qflux(n+u_ldown,l) *ne(n,ldown) & 279 + qflux(n+u_rup,l) *ne(n,rup) & 280 + qflux(n+u_left,l) *ne(n,left) & 281 + qflux(n+u_rdown,l) *ne(n,rdown) 282 283 newmass = mass(n,l) - dmass/Ai(n) 284 qi(n,l) = (qi(n,l)*mass(n,l) - dq/Ai(n) ) / newmass 285 IF(update_mass) mass(n,l) = newmass 286 END DO 256 dmass = hfluxt(ij+u_right,l)*ne(ij,right) & 257 + hfluxt(ij+u_lup,l) *ne(ij,lup) & 258 + hfluxt(ij+u_ldown,l)*ne(ij,ldown) & 259 + hfluxt(ij+u_rup,l) *ne(ij,rup) & 260 + hfluxt(ij+u_left,l) *ne(ij,left) & 261 + hfluxt(ij+u_rdown,l)*ne(ij,rdown) 262 263 dq = qflux(ij+u_right,l) *ne(ij,right) & 264 + qflux(ij+u_lup,l) *ne(ij,lup) & 265 + qflux(ij+u_ldown,l) *ne(ij,ldown) & 266 + qflux(ij+u_rup,l) *ne(ij,rup) & 267 + qflux(ij+u_left,l) *ne(ij,left) & 268 + qflux(ij+u_rdown,l) *ne(ij,rdown) 269 270 271 newmass = mass(ij,l) - dmass/Ai(ij) 272 qi(ij,l) = (qi(ij,l)*mass(ij,l) - dq/Ai(ij) ) / newmass 273 IF(update_mass) mass(ij,l) = newmass 274 287 275 END DO 288 276 END DO
Note: See TracChangeset
for help on using the changeset viewer.
