Changeset 174 for codes/icosagcm/trunk
- Timestamp:
- 10/07/13 17:48:24 (11 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 8 edited
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 -
codes/icosagcm/trunk/src/advect_tracer.f90
r151 r174 10 10 TYPE(t_field),POINTER :: f_one_over_sqrt_leng(:) 11 11 12 TYPE(t_message) :: req_u, req_ wfluxt, req_q, req_rhodz, req_gradq3d12 TYPE(t_message) :: req_u, req_cc, req_wfluxt, req_q, req_rhodz, req_gradq3d 13 13 14 14 REAL(rstd), PARAMETER :: pente_max=2.0 ! for vlz … … 81 81 first=.FALSE. 82 82 CALL init_message(f_u,req_e1_vect,req_u) 83 CALL init_message(f_cc,req_e1_scal,req_cc) 83 84 CALL init_message(f_wfluxt,req_i1,req_wfluxt) 84 85 CALL init_message(f_q,req_i1,req_q) … … 124 125 END DO 125 126 127 CALL send_message(f_cc,req_cc) 128 126 129 127 130 ! horizontal transport - split in two to place transfer of gradq3d … … 138 141 139 142 CALL send_message(f_gradq3d,req_gradq3d) 143 CALL wait_message(req_cc) 140 144 CALL wait_message(req_gradq3d) 141 145 … … 205 209 206 210 REAL(rstd) :: dzqmax, newmass, sigw, qq, w 207 INTEGER :: i,ij,l,j 211 INTEGER :: i,ij,l,j,ijb,ije 208 212 209 213 CALL trace_start("vlz") 214 215 ijb=((jj_begin-halo)-1)*iim+ii_begin-halo 216 ije = ((jj_end+halo)-1)*iim+ii_end+halo 210 217 211 218 ! finite difference of q 212 219 213 220 DO l=ll_beginp1,ll_end 214 DO j=jj_begin-halo,jj_end+halo 215 DO i=ii_begin-halo,ii_end+halo 216 ij=(j-1)*iim+i 217 dzqw(ij,l)=q(ij,l)-q(ij,l-1) 218 adzqw(ij,l)=abs(dzqw(ij,l)) 219 ENDDO 221 !$SIMD 222 DO ij=ijb,ije 223 dzqw(ij,l)=q(ij,l)-q(ij,l-1) 224 adzqw(ij,l)=abs(dzqw(ij,l)) 220 225 ENDDO 221 226 ENDDO … … 228 233 229 234 DO l=ll_beginp1,ll_endm1 230 DO j=jj_begin-halo,jj_end+halo 231 DO i=ii_begin-halo,ii_end+halo 232 ij=(j-1)*iim+i 233 IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN 234 dzq(ij,l) = 0.5*( dzqw(ij,l)+dzqw(ij,l+1) ) 235 dzqmax = pente_max * min( adzqw(ij,l),adzqw(ij,l+1) ) 236 dzq(ij,l) = sign( min(abs(dzq(ij,l)),dzqmax) , dzq(ij,l) ) ! NB : sign(a,b)=a*sign(b) 237 ELSE 238 dzq(ij,l)=0. 239 ENDIF 240 ENDDO 235 !$SIMD 236 DO ij=ijb,ije 237 IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN 238 dzq(ij,l) = 0.5*( dzqw(ij,l)+dzqw(ij,l+1) ) 239 dzqmax = pente_max * min( adzqw(ij,l),adzqw(ij,l+1) ) 240 dzq(ij,l) = sign( min(abs(dzq(ij,l)),dzqmax) , dzq(ij,l) ) ! NB : sign(a,b)=a*sign(b) 241 ELSE 242 dzq(ij,l)=0. 243 ENDIF 241 244 ENDDO 242 245 ENDDO … … 245 248 ! 0 slope in top and bottom layers 246 249 IF (omp_first) THEN 247 DO j=jj_begin-halo,jj_end+halo 248 DO i=ii_begin-halo,ii_end+halo 249 ij=(j-1)*iim+i 250 DO ij=ijb,ije 250 251 dzq(ij,1)=0. 251 ENDDO252 252 ENDDO 253 253 ENDIF 254 254 255 255 IF (omp_last) THEN 256 DO j=jj_begin-halo,jj_end+halo 257 DO i=ii_begin-halo,ii_end+halo 256 DO ij=ijb,ije 258 257 dzq(ij,llm)=0. 259 ENDDO260 258 ENDDO 261 259 ENDIF … … 267 265 ! then amount of q leaving level l/l+1 = wq = w * qq 268 266 DO l=ll_beginp1,ll_end 269 DO j=jj_begin-halo,jj_end+halo 270 DO i=ii_begin-halo,ii_end+halo 271 ij=(j-1)*iim+i 267 !$SIMD 268 DO ij=ijb,ije 272 269 w = fac*wfluxt(ij,l) 273 270 IF(w>0.) THEN ! upward transport, upwind side is at level l … … 279 276 ENDIF 280 277 wq(ij,l) = w*qq 281 ENDDO282 278 ENDDO 283 279 END DO 284 280 ! wq = 0 at top and bottom 285 281 IF (omp_first) THEN 286 DO j=jj_begin-halo,jj_end+halo 287 DO i=ii_begin-halo,ii_end+halo 288 ij=(j-1)*iim+i 282 DO ij=ijb,ije 289 283 wq(ij,1)=0. 290 ENDDO291 284 END DO 292 285 ENDIF 293 286 294 287 IF (omp_last) THEN 295 DO j=jj_begin-halo,jj_end+halo 296 DO i=ii_begin-halo,ii_end+halo 297 ij=(j-1)*iim+i 288 DO ij=ijb,ije 298 289 wq(ij,llm+1)=0. 299 ENDDO300 290 END DO 301 291 ENDIF … … 307 297 ! update q, mass is updated only after all q's have been updated 308 298 DO l=ll_begin,ll_end 309 DO j=jj_begin-halo,jj_end+halo 310 DO i=ii_begin-halo,ii_end+halo 311 ij=(j-1)*iim+i 299 !$SIMD 300 DO ij=ijb,ije 312 301 newmass = mass(ij,l) + fac*(wfluxt(ij,l)-wfluxt(ij,l+1)) 313 302 q(ij,l) = ( q(ij,l)*mass(ij,l) + wq(ij,l)-wq(ij,l+1) ) / newmass 314 303 IF(update_mass) mass(ij,l)=newmass 315 ENDDO316 304 ENDDO 317 305 END DO -
codes/icosagcm/trunk/src/caldyn_gcm.f90
r171 r174 97 97 wwuu=f_wwuu(ind) 98 98 99 DO j=jj_begin-1,jj_end+1 100 DO i=ii_begin-1,ii_end+1 101 ij=(j-1)*iim+i 102 ! lower BCs : geopot=phis, wflux=0, wwuu=0 103 geopot(ij,1) = phis(ij) 104 wflux(ij,1) = 0. 105 wwuu(ij+u_right,1)=0 106 wwuu(ij+u_lup,1)=0 107 wwuu(ij+u_ldown,1)=0 108 ! top BCs : wflux=0, wwuu=0 109 wflux(ij,llm+1) = 0. 110 wwuu(ij+u_right,llm+1)=0 111 wwuu(ij+u_lup,llm+1)=0 112 wwuu(ij+u_ldown,llm+1)=0 113 ENDDO 99 DO ij=ij_begin_ext,ij_end_ext 100 ! lower BCs : geopot=phis, wflux=0, wwuu=0 101 geopot(ij,1) = phis(ij) 102 wflux(ij,1) = 0. 103 wwuu(ij+u_right,1)=0 104 wwuu(ij+u_lup,1)=0 105 wwuu(ij+u_ldown,1)=0 106 ! top BCs : wflux=0, wwuu=0 107 wflux(ij,llm+1) = 0. 108 wwuu(ij+u_right,llm+1)=0 109 wwuu(ij+u_lup,llm+1)=0 110 wwuu(ij+u_ldown,llm+1)=0 114 111 ENDDO 115 112 END DO … … 303 300 ! REAL(rstd) :: qv(2*iim*jjm,llm) ! potential velocity 304 301 305 306 302 CALL trace_start("compute_pvort") 307 303 … … 316 312 DO l = ll_begin,ll_end 317 313 CALL test_message(req_u) 318 DO j=jj_begin-1,jj_end+1 319 DO i=ii_begin-1,ii_end+1 320 ij=(j-1)*iim+i 321 m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g 322 rhodz(ij,l) = m 323 theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 324 ENDDO 314 !DIR$ SIMD 315 DO ij=ij_begin_ext,ij_end_ext 316 m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g 317 rhodz(ij,l) = m 318 theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 325 319 ENDDO 326 320 ENDDO … … 328 322 DO l = ll_begin,ll_end 329 323 CALL test_message(req_u) 330 DO j=jj_begin-1,jj_end+1 331 DO i=ii_begin-1,ii_end+1 332 ij=(j-1)*iim+i 333 theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 334 ENDDO 335 ENDDO 324 !DIR$ SIMD 325 DO ij=ij_begin_ext,ij_end_ext 326 theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 327 ENDDO 336 328 ENDDO 337 329 END IF … … 341 333 !!! Compute shallow-water potential vorticity 342 334 DO l = ll_begin,ll_end 343 344 DO j=jj_begin-1,jj_end+1 345 DO i=ii_begin-1,ii_end+1 346 ij=(j-1)*iim+i 347 335 !DIR$ SIMD 336 DO ij=ij_begin_ext,ij_end_ext 348 337 etav= 1./Av(ij+z_up)*( ne_rup * u(ij+u_rup,l) * de(ij+u_rup) & 349 338 + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left) & … … 366 355 qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv 367 356 368 ENDDO369 357 ENDDO 370 358 371 DO j=jj_begin,jj_end 372 DO i=ii_begin,ii_end 373 ij=(j-1)*iim+i 374 qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l)) 375 qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l)) 376 qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l)) 377 END DO 359 !DIR$ SIMD 360 DO ij=ij_begin,ij_end 361 qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l)) 362 qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l)) 363 qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l)) 378 364 END DO 379 365 … … 407 393 DO l = 1,llm 408 394 !$OMP DO SCHEDULE(STATIC) 409 DO j=jj_begin-1,jj_end+1 410 DO i=ii_begin-1,ii_end+1 411 ij=(j-1)*iim+i 412 p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later 413 ! p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j)) 414 exner_ik = cpp * (p_ik/preff) ** kappa 415 pk(ij,l) = exner_ik 416 ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 417 geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 418 ENDDO 395 !DIR$ SIMD 396 DO ij=ij_begin_ext,ij_end_ext 397 p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later 398 ! p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j)) 399 exner_ik = cpp * (p_ik/preff) ** kappa 400 pk(ij,l) = exner_ik 401 ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 402 geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 419 403 ENDDO 420 404 ENDDO … … 427 411 428 412 ! uppermost layer 429 DO j=jj_begin-1,jj_end+1 430 DO i=ii_begin-1,ii_end+1 431 ij=(j-1)*iim+i 432 pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm) 433 END DO 413 !DIR$ SIMD 414 DO ij=ij_begin_ext,ij_end_ext 415 pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm) 434 416 END DO 435 417 ! other layers 436 418 DO l = llm-1, 1, -1 437 419 !$OMP DO SCHEDULE(STATIC) 438 DO j=jj_begin-1,jj_end+1 439 DO i=ii_begin-1,ii_end+1 440 ij=(j-1)*iim+i 441 pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1)) 442 END DO 420 !DIR$ SIMD 421 DO ij=ij_begin_ext,ij_end_ext 422 pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1)) 443 423 END DO 444 424 END DO 445 425 ! surface pressure (for diagnostics) 446 DO j=jj_begin-1,jj_end+1 447 DO i=ii_begin-1,ii_end+1 448 ij=(j-1)*iim+i 449 ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1) 450 END DO 426 DO ij=ij_begin_ext,ij_end_ext 427 ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1) 451 428 END DO 452 429 … … 454 431 DO l = 1,llm 455 432 !$OMP DO SCHEDULE(STATIC) 456 DO j=jj_begin-1,jj_end+1 457 DO i=ii_begin-1,ii_end+1 458 ij=(j-1)*iim+i 459 geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l) 460 ENDDO 433 !DIR$ SIMD 434 DO ij=ij_begin_ext,ij_end_ext 435 geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l) 461 436 ENDDO 462 437 ENDDO … … 464 439 DO l = 1,llm 465 440 !$OMP DO SCHEDULE(STATIC) 466 DO j=jj_begin-1,jj_end+1 467 DO i=ii_begin-1,ii_end+1 468 ij=(j-1)*iim+i 469 p_ik = pk(ij,l) 470 exner_ik = cpp * (p_ik/preff) ** kappa 471 geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 472 pk(ij,l) = exner_ik 473 ENDDO 441 !DIR$ SIMD 442 DO ij=ij_begin_ext,ij_end_ext 443 p_ik = pk(ij,l) 444 exner_ik = cpp * (p_ik/preff) ** kappa 445 geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 446 pk(ij,l) = exner_ik 474 447 ENDDO 475 448 ENDDO … … 514 487 !!! Compute mass and theta fluxes 515 488 IF (caldyn_conserv==energy) CALL test_message(req_qu) 516 DO j=jj_begin-1,jj_end+1 517 DO i=ii_begin-1,ii_end+1 518 ij=(j-1)*iim+i 489 !DIR$ SIMD 490 DO ij=ij_begin_ext,ij_end_ext 519 491 hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) 520 492 hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) … … 524 496 Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l) 525 497 Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*hflux(ij+u_ldown,l) 526 ENDDO527 498 ENDDO 528 499 529 500 !!! compute horizontal divergence of fluxes 530 DO j=jj_begin,jj_end 531 DO i=ii_begin,ii_end 532 ij=(j-1)*iim+i 501 !DIR$ SIMD 502 DO ij=ij_begin,ij_end 533 503 ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 534 504 convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l) + & … … 547 517 ne_ldown*Ftheta(ij+u_ldown,l) + & 548 518 ne_rdown*Ftheta(ij+u_rdown,l)) 549 ENDDO550 519 ENDDO 551 520 … … 560 529 561 530 DO l=ll_begin,ll_end 562 DO j=jj_begin,jj_end 563 DO i=ii_begin,ii_end 564 ij=(j-1)*iim+i 565 531 !DIR$ SIMD 532 DO ij=ij_begin,ij_end 533 566 534 uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+ & 567 535 wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+ & … … 601 569 du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown) 602 570 603 ENDDO604 571 ENDDO 605 572 ENDDO … … 608 575 609 576 DO l=ll_begin,ll_end 610 DO j=jj_begin,jj_end 611 DO i=ii_begin,ii_end 612 ij=(j-1)*iim+i 577 !DIR$ SIMD 578 DO ij=ij_begin,ij_end 613 579 614 580 uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+ & … … 649 615 du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown) 650 616 651 ENDDO652 617 ENDDO 653 618 ENDDO … … 661 626 662 627 DO l=ll_begin,ll_end 663 DO j=jj_begin,jj_end 664 DO i=ii_begin,ii_end 665 ij=(j-1)*iim+i 628 !DIR$ SIMD 629 DO ij=ij_begin,ij_end 666 630 667 631 berni(ij,l) = pk(ij,l) + & … … 673 637 le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 674 638 pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) 675 ENDDO676 639 ENDDO 677 640 ENDDO … … 680 643 681 644 DO l=ll_begin,ll_end 682 DO j=jj_begin,jj_end 683 DO i=ii_begin,ii_end 684 ij=(j-1)*iim+i 645 !DIR$ SIMD 646 DO ij=ij_begin,ij_end 685 647 686 648 berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & … … 691 653 le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 + & 692 654 le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 693 694 ENDDO695 655 ENDDO 696 656 ENDDO … … 700 660 !!! Add gradients of Bernoulli and Exner functions to du 701 661 DO l=ll_begin,ll_end 702 DO j=jj_begin,jj_end 703 DO i=ii_begin,ii_end 704 ij=(j-1)*iim+i 662 !DIR$ SIMD 663 DO ij=ij_begin,ij_end 705 664 706 665 du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * ( & … … 720 679 + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) ) 721 680 722 ENDDO723 681 ENDDO 724 682 ENDDO … … 750 708 REAL(rstd) :: p_ik, exner_ik 751 709 710 752 711 CALL trace_start("compute_caldyn_vert") 753 712 … … 757 716 IF (caldyn_conserv==energy) CALL test_message(req_qu) 758 717 !$OMP DO SCHEDULE(STATIC) 759 DO j=jj_begin,jj_end 760 DO i=ii_begin,ii_end 761 ij=(j-1)*iim+i 718 !DIR$ SIMD 719 DO ij=ij_begin,ij_end 762 720 convm(ij,l) = convm(ij,l) + convm(ij,l+1) 763 ENDDO764 721 ENDDO 765 722 ENDDO … … 770 727 ! compute dps 771 728 IF (omp_first) THEN 772 DO j=jj_begin,jj_end 773 DO i=ii_begin,ii_end 774 ij=(j-1)*iim+i 729 !DIR$ SIMD 730 DO ij=ij_begin,ij_end 775 731 ! dps/dt = -int(div flux)dz 776 732 dps(ij) = convm(ij,1) * g 777 ENDDO778 733 ENDDO 779 734 ENDIF … … 782 737 DO l=ll_beginp1,ll_end 783 738 IF (caldyn_conserv==energy) CALL test_message(req_qu) 784 DO j=jj_begin,jj_end 785 DO i=ii_begin,ii_end 786 ij=(j-1)*iim+i 739 !DIR$ SIMD 740 DO ij=ij_begin,ij_end 787 741 ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt 788 742 ! => w>0 for upward transport 789 743 wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l ) 790 ENDDO791 744 ENDDO 792 745 ENDDO 793 746 794 747 DO l=ll_begin,ll_endm1 795 DO j=jj_begin,jj_end 796 DO i=ii_begin,ii_end 797 ij=(j-1)*iim+i 748 !DIR$ SIMD 749 DO ij=ij_begin,ij_end 798 750 dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l ) - 0.5 * ( wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1))) 799 ENDDO800 751 ENDDO 801 752 ENDDO 802 753 803 754 DO l=ll_beginp1,ll_end 804 DO j=jj_begin,jj_end 805 DO i=ii_begin,ii_end 806 ij=(j-1)*iim+i 755 !DIR$ SIMD 756 DO ij=ij_begin,ij_end 807 757 dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l ) + 0.5 * ( wflux(ij,l ) * (theta(ij,l-1) + theta(ij,l) ) ) 808 ENDDO809 758 ENDDO 810 759 ENDDO … … 812 761 ! Compute vertical transport 813 762 DO l=ll_beginp1,ll_end 814 DO j=jj_begin,jj_end 815 DO i=ii_begin,ii_end 816 ij=(j-1)*iim+i 763 !DIR$ SIMD 764 DO ij=ij_begin,ij_end 817 765 wwuu(ij+u_right,l) = 0.5*( wflux(ij,l) + wflux(ij+t_right,l)) * (u(ij+u_right,l) - u(ij+u_right,l-1)) 818 766 wwuu(ij+u_lup,l) = 0.5* ( wflux(ij,l) + wflux(ij+t_lup,l)) * (u(ij+u_lup,l) - u(ij+u_lup,l-1)) 819 767 wwuu(ij+u_ldown,l) = 0.5*( wflux(ij,l) + wflux(ij+t_ldown,l)) * (u(ij+u_ldown,l) - u(ij+u_ldown,l-1)) 820 ENDDO821 768 ENDDO 822 769 ENDDO … … 827 774 ! Add vertical transport to du 828 775 DO l=ll_begin,ll_end 829 DO j=jj_begin,jj_end 830 DO i=ii_begin,ii_end 831 ij=(j-1)*iim+i 776 !DIR$ SIMD 777 DO ij=ij_begin,ij_end 832 778 du(ij+u_right, l ) = du(ij+u_right,l) - (wwuu(ij+u_right,l+1)+ wwuu(ij+u_right,l)) / (rhodz(ij,l)+rhodz(ij+t_right,l)) 833 779 du(ij+u_lup, l ) = du(ij+u_lup,l) - (wwuu(ij+u_lup,l+1) + wwuu(ij+u_lup,l)) / (rhodz(ij,l)+rhodz(ij+t_lup,l)) 834 780 du(ij+u_ldown, l ) = du(ij+u_ldown,l) - (wwuu(ij+u_ldown,l+1)+ wwuu(ij+u_ldown,l)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l)) 835 ENDDO836 781 ENDDO 837 782 ENDDO -
codes/icosagcm/trunk/src/dimensions.f90
r151 r174 9 9 INTEGER :: ii_nb 10 10 INTEGER :: jj_nb 11 11 INTEGER :: ij_begin 12 INTEGER :: ij_end 13 INTEGER :: ij_begin_ext 14 INTEGER :: ij_end_ext 15 12 16 INTEGER :: t_right 13 17 INTEGER :: t_rup … … 56 60 ii_nb=d%ii_nb 57 61 jj_nb=d%jj_nb 58 62 63 ij_begin=(jj_begin-1)*iim+ii_begin 64 ij_end=(jj_end-1)*iim+ii_end 65 ij_begin_ext=((jj_begin-1)-1)*iim+ii_begin-1 66 ij_end_ext = ((jj_end+1)-1)*iim+ii_end+1 67 59 68 t_right=d%t_right 60 69 t_rup=d%t_rup -
codes/icosagcm/trunk/src/dissip_gcm.f90
r167 r174 77 77 CHARACTER(len=255) :: rayleigh_friction_key 78 78 REAL(rstd) :: mintau 79 INTEGER :: seed_size 80 INTEGER,ALLOCATABLE :: seed(:) 79 81 80 82 81 INTEGER :: i,j, n,ind,it,iter83 INTEGER :: i,j,ij,ind,it,iter 82 84 83 85 rayleigh_friction_key='none' … … 142 144 cdivgrad=-1 143 145 cgradrot=-1 144 145 CALL RANDOM_SEED()146 146 147 147 DO ind=1,ndomain … … 152 152 DO j=jj_begin,jj_end 153 153 DO i=ii_begin,ii_end 154 n=(j-1)*iim+i154 ij=(j-1)*iim+i 155 155 CALL RANDOM_NUMBER(r) 156 u( n+u_right)=r-0.5156 u(ij+u_right)=r-0.5 157 157 CALL RANDOM_NUMBER(r) 158 u( n+u_lup)=r-0.5158 u(ij+u_lup)=r-0.5 159 159 CALL RANDOM_NUMBER(r) 160 u(n+u_ldown)=r-0.5 161 ENDDO 162 ENDDO 163 160 u(ij+u_ldown)=r-0.5 161 ENDDO 162 ENDDO 164 163 ENDDO 165 164 … … 191 190 DO j=jj_begin,jj_end 192 191 DO i=ii_begin,ii_end 193 n=(j-1)*iim+i194 if (le( n+u_right)>1e-100) dumax=MAX(dumax,ABS(du(n+u_right)))195 if (le( n+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(n+u_lup)))196 if (le( n+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(n+u_ldown)))192 ij=(j-1)*iim+i 193 if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right))) 194 if (le(ij+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_lup))) 195 if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown))) 197 196 ENDDO 198 197 ENDDO … … 228 227 u=f_u(ind) 229 228 230 DO j=jj_begin,jj_end229 DO j=jj_begin,jj_end 231 230 DO i=ii_begin,ii_end 232 n=(j-1)*iim+i231 ij=(j-1)*iim+i 233 232 CALL RANDOM_NUMBER(r) 234 u( n+u_right)=r-0.5233 u(ij+u_right)=r-0.5 235 234 CALL RANDOM_NUMBER(r) 236 u( n+u_lup)=r-0.5235 u(ij+u_lup)=r-0.5 237 236 CALL RANDOM_NUMBER(r) 238 u(n+u_ldown)=r-0.5 239 ENDDO 240 ENDDO 241 237 u(ij+u_ldown)=r-0.5 238 ENDDO 239 ENDDO 242 240 ENDDO 243 241 … … 268 266 DO j=jj_begin,jj_end 269 267 DO i=ii_begin,ii_end 270 n=(j-1)*iim+i271 if (le( n+u_right)>1e-100) dumax=MAX(dumax,ABS(du(n+u_right)))272 if (le( n+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(n+u_lup)))273 if (le( n+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(n+u_ldown)))268 ij=(j-1)*iim+i 269 if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right))) 270 if (le(ij+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_lup))) 271 if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown))) 274 272 ENDDO 275 273 ENDDO … … 304 302 theta=f_theta(ind) 305 303 306 DO j=jj_begin,jj_end 307 DO i=ii_begin,ii_end 308 n=(j-1)*iim+i 309 CALL RANDOM_NUMBER(r) 310 theta(n)=r-0.5 311 ENDDO 312 ENDDO 313 304 DO j=jj_begin,jj_end 305 DO i=ii_begin,ii_end 306 ij=(j-1)*iim+i 307 CALL RANDOM_NUMBER(r) 308 theta(ij)=r-0.5 309 ENDDO 310 ENDDO 314 311 ENDDO 315 312 … … 339 336 DO j=jj_begin,jj_end 340 337 DO i=ii_begin,ii_end 341 n=(j-1)*iim+i342 dthetamax=MAX(dthetamax,ABS(dtheta( n)))338 ij=(j-1)*iim+i 339 dthetamax=MAX(dthetamax,ABS(dtheta(ij))) 343 340 ENDDO 344 341 ENDDO 345 342 ENDDO 343 346 344 IF (using_mpi) THEN 347 345 CALL MPI_ALLREDUCE(dthetamax,dthetamax1,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 348 346 dthetamax=dthetamax1 349 347 ENDIF 348 350 349 IF (is_mpi_root) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax 351 350 … … 417 416 418 417 INTEGER :: ind, shear 419 INTEGER :: l,i ,j,n418 INTEGER :: l,ij 420 419 421 420 !$OMP BARRIER … … 444 443 445 444 DO l=ll_begin,ll_end 446 DO j=jj_begin,jj_end 447 DO i=ii_begin,ii_end 448 n=(j-1)*iim+i 449 450 due(n+u_right,l) = -0.5*( due_diss1(n+u_right,l)/tau_graddiv(l) + due_diss2(n+u_right,l)/tau_gradrot(l))*itau_dissip 451 due(n+u_lup,l) = -0.5*( due_diss1(n+u_lup,l) /tau_graddiv(l) + due_diss2(n+u_lup,l) /tau_gradrot(l))*itau_dissip 452 due(n+u_ldown,l) = -0.5*( due_diss1(n+u_ldown,l)/tau_graddiv(l) + due_diss2(n+u_ldown,l)/tau_gradrot(l))*itau_dissip 453 454 dtheta_rhodz(n,l) = -0.5*dtheta_rhodz_diss(n,l)/tau_divgrad(l)*itau_dissip 455 ENDDO 445 !$SIMD 446 DO ij=ij_begin,ij_end 447 448 due(ij+u_right,l) = -0.5*( due_diss1(ij+u_right,l)/tau_graddiv(l) + due_diss2(ij+u_right,l)/tau_gradrot(l))*itau_dissip 449 due(ij+u_lup,l) = -0.5*( due_diss1(ij+u_lup,l) /tau_graddiv(l) + due_diss2(ij+u_lup,l) /tau_gradrot(l))*itau_dissip 450 due(ij+u_ldown,l) = -0.5*( due_diss1(ij+u_ldown,l)/tau_graddiv(l) + due_diss2(ij+u_ldown,l)/tau_gradrot(l))*itau_dissip 451 452 dtheta_rhodz(ij,l) = -0.5*dtheta_rhodz_diss(ij,l)/tau_divgrad(l)*itau_dissip 456 453 ENDDO 457 454 ENDDO … … 490 487 LOGICAL :: hybrid_eta 491 488 INTEGER :: shift_u, shift_t, zcoords, nn 492 z = (phi( n,l)+phi(n+shift_t,l))/(2.*g)489 z = (phi(ij,l)+phi(ij+shift_t,l))/(2.*g) 493 490 IF(z>zh) THEN ! relax only in the sponge layer z>zh 494 491 495 nn = n+shift_u492 nn = ij+shift_u 496 493 CALL xyz2lonlat(xyz_e(nn,:),lon,lat) 497 494 zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm … … 522 519 REAL(rstd),POINTER :: due(:,:) 523 520 INTEGER :: ind 524 INTEGER :: it, i,j,l,ij521 INTEGER :: it,l,ij 525 522 526 523 CALL trace_start("gradiv") … … 532 529 due=f_due(ind) 533 530 DO l = ll_begin, ll_end 534 DO j=jj_begin,jj_end 535 DO i=ii_begin,ii_end 536 ij=(j-1)*iim+i 531 !$SIMD 532 DO ij=ij_begin,ij_end 537 533 due(ij+u_right,l)=ue(ij+u_right,l) 538 534 due(ij+u_lup,l)=ue(ij+u_lup,l) 539 535 due(ij+u_ldown,l)=ue(ij+u_ldown,l) 540 ENDDO541 536 ENDDO 542 537 ENDDO … … 571 566 REAL(rstd),POINTER :: due(:,:) 572 567 INTEGER :: ind 573 INTEGER :: it, i,j,l,ij568 INTEGER :: it,l,ij 574 569 575 570 CALL trace_start("gradrot") … … 581 576 due=f_due(ind) 582 577 DO l = ll_begin, ll_end 583 DO j=jj_begin,jj_end 584 DO i=ii_begin,ii_end 585 ij=(j-1)*iim+i 578 !$SIMD 579 DO ij=ij_begin,ij_end 586 580 due(ij+u_right,l)=ue(ij+u_right,l) 587 581 due(ij+u_lup,l)=ue(ij+u_lup,l) 588 582 due(ij+u_ldown,l)=ue(ij+u_ldown,l) 589 ENDDO590 583 ENDDO 591 584 ENDDO … … 663 656 664 657 INTEGER :: ind 665 INTEGER :: it, i,j,l,ij658 INTEGER :: it,l,ij 666 659 667 660 CALL trace_start("divgrad") … … 674 667 dtheta_rhodz=f_dtheta_rhodz(ind) 675 668 DO l = ll_begin, ll_end 676 DO j=jj_begin,jj_end 677 DO i=ii_begin,ii_end 678 ij=(j-1)*iim+i 669 !$SIMD 670 DO ij=ij_begin,ij_end 679 671 dtheta_rhodz(ij,l) = theta_rhodz(ij,l) / mass(ij,l) 680 ENDDO681 672 ENDDO 682 673 ENDDO … … 702 693 703 694 DO l = ll_begin, ll_end 704 DO j=jj_begin,jj_end 705 DO i=ii_begin,ii_end 706 ij=(j-1)*iim+i 695 !$SIMD 696 DO ij=ij_begin,ij_end 707 697 dtheta_rhodz(ij,l) = dtheta_rhodz(ij,l) * mass(ij,l) 708 ENDDO709 698 ENDDO 710 699 ENDDO … … 725 714 REAL(rstd) :: divu_i(iim*jjm,llb:lle) 726 715 727 INTEGER :: i ,j,n,l716 INTEGER :: ij,l 728 717 729 718 DO l=llb,lle 730 DO j=jj_begin,jj_end 731 DO i=ii_begin,ii_end 732 n=(j-1)*iim+i 733 divu_i(n,l)=1./Ai(n)*(ne(n,right)*ue(n+u_right,l)*le(n+u_right) + & 734 ne(n,rup)*ue(n+u_rup,l)*le(n+u_rup) + & 735 ne(n,lup)*ue(n+u_lup,l)*le(n+u_lup) + & 736 ne(n,left)*ue(n+u_left,l)*le(n+u_left) + & 737 ne(n,ldown)*ue(n+u_ldown,l)*le(n+u_ldown) + & 738 ne(n,rdown)*ue(n+u_rdown,l)*le(n+u_rdown)) 739 ENDDO 719 !$SIMD 720 DO ij=ij_begin,ij_end 721 divu_i(ij,l)=1./Ai(ij)*(ne(ij,right)*ue(ij+u_right,l)*le(ij+u_right) + & 722 ne(ij,rup)*ue(ij+u_rup,l)*le(ij+u_rup) + & 723 ne(ij,lup)*ue(ij+u_lup,l)*le(ij+u_lup) + & 724 ne(ij,left)*ue(ij+u_left,l)*le(ij+u_left) + & 725 ne(ij,ldown)*ue(ij+u_ldown,l)*le(ij+u_ldown) + & 726 ne(ij,rdown)*ue(ij+u_rdown,l)*le(ij+u_rdown)) 740 727 ENDDO 741 728 ENDDO 742 729 743 730 DO l=llb,lle 744 DO j=jj_begin,jj_end 745 DO i=ii_begin,ii_end731 !$SIMD 732 DO ij=ij_begin,ij_end 746 733 747 n=(j-1)*iim+i 748 749 gradivu_e(n+u_right,l)=-1/de(n+u_right)*(ne(n,right)*divu_i(n,l)+ ne(n+t_right,left)*divu_i(n+t_right,l) ) 750 751 gradivu_e(n+u_lup,l)=-1/de(n+u_lup)*(ne(n,lup)*divu_i(n,l)+ ne(n+t_lup,rdown)*divu_i(n+t_lup,l)) 752 753 gradivu_e(n+u_ldown,l)=-1/de(n+u_ldown)*(ne(n,ldown)*divu_i(n,l)+ne(n+t_ldown,rup)*divu_i(n+t_ldown,l) ) 754 755 ENDDO 734 gradivu_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*divu_i(ij,l)+ ne(ij+t_right,left)*divu_i(ij+t_right,l) ) 735 736 gradivu_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*divu_i(ij,l)+ ne(ij+t_lup,rdown)*divu_i(ij+t_lup,l)) 737 738 gradivu_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*divu_i(ij,l)+ne(ij+t_ldown,rup)*divu_i(ij+t_ldown,l) ) 739 756 740 ENDDO 757 741 ENDDO 758 742 759 743 DO l=llb,lle 760 DO j=jj_begin,jj_end 761 DO i=ii_begin,ii_end 762 n=(j-1)*iim+i 763 gradivu_e(n+u_right,l)=-gradivu_e(n+u_right,l)*cgraddiv 764 gradivu_e(n+u_lup,l)=-gradivu_e(n+u_lup,l)*cgraddiv 765 gradivu_e(n+u_ldown,l)=-gradivu_e(n+u_ldown,l)*cgraddiv 766 ENDDO 744 !$SIMD 745 DO ij=ij_begin,ij_end 746 gradivu_e(ij+u_right,l)=-gradivu_e(ij+u_right,l)*cgraddiv 747 gradivu_e(ij+u_lup,l)=-gradivu_e(ij+u_lup,l)*cgraddiv 748 gradivu_e(ij+u_ldown,l)=-gradivu_e(ij+u_ldown,l)*cgraddiv 767 749 ENDDO 768 750 ENDDO … … 780 762 REAL(rstd) :: grad_e(3*iim*jjm,llb:lle) 781 763 782 INTEGER :: i ,j,n,l764 INTEGER :: ij,l 783 765 784 766 785 767 DO l=llb,lle 786 DO j=jj_begin-1,jj_end+1 787 DO i=ii_begin-1,ii_end+1 788 789 n=(j-1)*iim+i 768 !$SIMD 769 DO ij=ij_begin_ext,ij_end_ext 790 770 791 grad_e(n+u_right,l)=-1/de(n+u_right)*(ne(n,right)*theta(n,l)+ ne(n+t_right,left)*theta(n+t_right,l) ) 792 793 grad_e(n+u_lup,l)=-1/de(n+u_lup)*(ne(n,lup)*theta(n,l)+ ne(n+t_lup,rdown)*theta(n+t_lup,l )) 794 795 grad_e(n+u_ldown,l)=-1/de(n+u_ldown)*(ne(n,ldown)*theta(n,l)+ne(n+t_ldown,rup)*theta(n+t_ldown,l) ) 796 797 ENDDO 771 grad_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*theta(ij,l)+ ne(ij+t_right,left)*theta(ij+t_right,l) ) 772 773 grad_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*theta(ij,l)+ ne(ij+t_lup,rdown)*theta(ij+t_lup,l )) 774 775 grad_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*theta(ij,l)+ne(ij+t_ldown,rup)*theta(ij+t_ldown,l) ) 776 798 777 ENDDO 799 778 ENDDO … … 801 780 802 781 DO l=llb,lle 803 DO j=jj_begin,jj_end 804 DO i=ii_begin,ii_end 805 n=(j-1)*iim+i 806 divgrad_i(n,l)=1./Ai(n)*(ne(n,right)*grad_e(n+u_right,l)*le(n+u_right) + & 807 ne(n,rup)*grad_e(n+u_rup,l)*le(n+u_rup) + & 808 ne(n,lup)*grad_e(n+u_lup,l)*le(n+u_lup) + & 809 ne(n,left)*grad_e(n+u_left,l)*le(n+u_left) + & 810 ne(n,ldown)*grad_e(n+u_ldown,l)*le(n+u_ldown) + & 811 ne(n,rdown)*grad_e(n+u_rdown,l)*le(n+u_rdown)) 812 ENDDO 782 !$SIMD 783 DO ij=jj_begin,ij_end 784 785 divgrad_i(ij,l)=1./Ai(ij)*(ne(ij,right)*grad_e(ij+u_right,l)*le(ij+u_right) + & 786 ne(ij,rup)*grad_e(ij+u_rup,l)*le(ij+u_rup) + & 787 ne(ij,lup)*grad_e(ij+u_lup,l)*le(ij+u_lup) + & 788 ne(ij,left)*grad_e(ij+u_left,l)*le(ij+u_left) + & 789 ne(ij,ldown)*grad_e(ij+u_ldown,l)*le(ij+u_ldown) + & 790 ne(ij,rdown)*grad_e(ij+u_rdown,l)*le(ij+u_rdown)) 813 791 ENDDO 814 792 ENDDO 815 793 816 794 DO l=llb,lle 817 DO j=jj_begin,jj_end 818 DO i=ii_begin,ii_end 819 n=(j-1)*iim+i 820 divgrad_i(n,l)=-divgrad_i(n,l)*cdivgrad 821 ENDDO 795 DO ij=ij_begin,ij_end 796 divgrad_i(ij,l)=-divgrad_i(ij,l)*cdivgrad 822 797 ENDDO 823 798 ENDDO … … 835 810 REAL(rstd) :: rot_v(2*iim*jjm,llb:lle) 836 811 837 INTEGER :: i ,j,n,l812 INTEGER :: ij,l 838 813 839 814 DO l=llb,lle 840 DO j=jj_begin-1,jj_end+1 841 DO i=ii_begin-1,ii_end+1 842 n=(j-1)*iim+i 843 844 rot_v(n+z_up,l)= 1./Av(n+z_up)*( ne(n,rup)*ue(n+u_rup,l)*de(n+u_rup) & 845 + ne(n+t_rup,left)*ue(n+t_rup+u_left,l)*de(n+t_rup+u_left) & 846 - ne(n,lup)*ue(n+u_lup,l)*de(n+u_lup) ) 815 !$SIMD 816 DO ij=ij_begin_ext,ij_end_ext 817 818 rot_v(ij+z_up,l)= 1./Av(ij+z_up)*( ne(ij,rup)*ue(ij+u_rup,l)*de(ij+u_rup) & 819 + ne(ij+t_rup,left)*ue(ij+t_rup+u_left,l)*de(ij+t_rup+u_left) & 820 - ne(ij,lup)*ue(ij+u_lup,l)*de(ij+u_lup) ) 847 821 848 rot_v(n+z_down,l) = 1./Av(n+z_down)*( ne(n,ldown)*ue(n+u_ldown,l)*de(n+u_ldown) & 849 + ne(n+t_ldown,right)*ue(n+t_ldown+u_right,l)*de(n+t_ldown+u_right) & 850 - ne(n,rdown)*ue(n+u_rdown,l)*de(n+u_rdown) ) 851 852 ENDDO 822 rot_v(ij+z_down,l) = 1./Av(ij+z_down)*( ne(ij,ldown)*ue(ij+u_ldown,l)*de(ij+u_ldown) & 823 + ne(ij+t_ldown,right)*ue(ij+t_ldown+u_right,l)*de(ij+t_ldown+u_right) & 824 - ne(ij,rdown)*ue(ij+u_rdown,l)*de(ij+u_rdown) ) 825 853 826 ENDDO 854 827 ENDDO 855 828 856 829 DO l=llb,lle 857 DO j=jj_begin,jj_end 858 DO i=ii_begin,ii_end 859 n=(j-1)*iim+i 860 861 gradrot_e(n+u_right,l)=1/le(n+u_right)*ne(n,right)*(rot_v(n+z_rdown,l)-rot_v(n+z_rup,l)) 830 !$SIMD 831 DO ij=ij_begin,ij_end 832 833 gradrot_e(ij+u_right,l)=1/le(ij+u_right)*ne(ij,right)*(rot_v(ij+z_rdown,l)-rot_v(ij+z_rup,l)) 862 834 863 gradrot_e(n+u_lup,l)=1/le(n+u_lup)*ne(n,lup)*(rot_v(n+z_up,l)-rot_v(n+z_lup,l)) 864 865 gradrot_e(n+u_ldown,l)=1/le(n+u_ldown)*ne(n,ldown)*(rot_v(n+z_ldown,l)-rot_v(n+z_down,l)) 866 867 ENDDO 835 gradrot_e(ij+u_lup,l)=1/le(ij+u_lup)*ne(ij,lup)*(rot_v(ij+z_up,l)-rot_v(ij+z_lup,l)) 836 837 gradrot_e(ij+u_ldown,l)=1/le(ij+u_ldown)*ne(ij,ldown)*(rot_v(ij+z_ldown,l)-rot_v(ij+z_down,l)) 838 868 839 ENDDO 869 840 ENDDO 870 841 871 842 DO l=llb,lle 872 DO j=jj_begin,jj_end 873 DO i=ii_begin,ii_end 874 n=(j-1)*iim+i 875 876 gradrot_e(n+u_right,l)=-gradrot_e(n+u_right,l)*cgradrot 877 gradrot_e(n+u_lup,l)=-gradrot_e(n+u_lup,l)*cgradrot 878 gradrot_e(n+u_ldown,l)=-gradrot_e(n+u_ldown,l)*cgradrot 879 880 ENDDO 843 !$SIMD 844 DO ij=ij_begin,ij_end 845 846 gradrot_e(ij+u_right,l)=-gradrot_e(ij+u_right,l)*cgradrot 847 gradrot_e(ij+u_lup,l)=-gradrot_e(ij+u_lup,l)*cgradrot 848 gradrot_e(ij+u_ldown,l)=-gradrot_e(ij+u_ldown,l)*cgradrot 849 881 850 ENDDO 882 851 ENDDO -
codes/icosagcm/trunk/src/domain.f90
r151 r174 72 72 USE grid_param 73 73 USE mpipara 74 USE IOIPSL 74 75 IMPLICIT NONE 75 76 INTEGER :: ind,nf,ni,nj,i,j 76 77 INTEGER :: quotient, rest 78 INTEGER :: halo_i,halo_j 79 77 80 TYPE(t_domain),POINTER :: d 78 81 … … 81 84 ALLOCATE(domglo_rank(ndomain_glo)) 82 85 ALLOCATE(domglo_loc_ind(ndomain_glo)) 86 87 halo_i=0 88 CALL getin("halo_i",halo_i) 89 halo_j=1 90 CALL getin("halo_j",halo_j) 83 91 84 92 ind=0 … … 122 130 123 131 124 d%iim=d%ii_nb+2*halo 125 d%jjm=d%jj_nb+2*halo 126 d%ii_begin=halo+1 127 d%jj_begin=halo+1 132 d%iim=d%ii_nb+2*halo + halo_i*2 133 d%jjm=d%jj_nb+2*halo + halo_j*2 134 d%ii_begin=halo+1 + halo_i 135 d%jj_begin=halo+1 + halo_j 128 136 d%ii_end=d%ii_begin+d%ii_nb-1 129 137 d%jj_end=d%jj_begin+d%jj_nb-1 -
codes/icosagcm/trunk/src/domain_param.f90
r95 r174 3 3 INTEGER :: nsplit_i 4 4 INTEGER :: nsplit_j 5 INTEGER :: halo= 35 INTEGER :: halo=1 6 6 7 7 INTEGER, PARAMETER :: default_nsplit_i=1 -
codes/icosagcm/trunk/src/timeloop_gcm.f90
r171 r174 336 336 ps=f_ps(ind) ; dps=f_dps(ind) ; 337 337 IF (omp_first) THEN 338 DO j=jj_begin,jj_end 339 DO i=ii_begin,ii_end 340 ij=(j-1)*iim+i 341 ps(ij)=ps(ij)+dt*dps(ij) 342 ENDDO 338 !$SIMD 339 DO ij=ij_begin,ij_end 340 ps(ij)=ps(ij)+dt*dps(ij) 343 341 ENDDO 344 342 ENDIF … … 346 344 mass=f_mass(ind) ; dmass=f_dmass(ind) ; 347 345 DO l=1,llm 348 DO j=jj_begin,jj_end 349 DO i=ii_begin,ii_end 350 ij=(j-1)*iim+i 351 mass(ij,l)=mass(ij,l)+dt*dmass(ij,l) 352 ENDDO 346 !$SIMD 347 DO ij=ij_begin,ij_end 348 mass(ij,l)=mass(ij,l)+dt*dmass(ij,l) 353 349 ENDDO 354 350 END DO … … 364 360 365 361 DO l=ll_begin,ll_end 366 DO j=jj_begin,jj_end 367 DO i=ii_begin,ii_end 368 ij=(j-1)*iim+i 362 !$SIMD 363 DO ij=ij_begin,ij_end 369 364 u(ij+u_right,l)=u(ij+u_right,l)+dt*du(ij+u_right,l) 370 365 u(ij+u_lup,l)=u(ij+u_lup,l)+dt*du(ij+u_lup,l) 371 366 u(ij+u_ldown,l)=u(ij+u_ldown,l)+dt*du(ij+u_ldown,l) 372 367 theta_rhodz(ij,l)=theta_rhodz(ij,l)+dt*dtheta_rhodz(ij,l) 373 ENDDO374 368 ENDDO 375 369 ENDDO … … 401 395 402 396 IF (stage==1) THEN ! first stage : save model state in XXm1 403 DO j=jj_begin,jj_end 404 DO i=ii_begin,ii_end 405 ij=(j-1)*iim+i 406 psm1(ij)=ps(ij) 407 ENDDO 408 ENDDO 397 !$SIMD 398 DO ij=ij_begin,ij_end 399 psm1(ij)=ps(ij) 400 ENDDO 409 401 ENDIF 410 402 411 403 ! updates are of the form x1 := x0 + tau*f(x1) 412 DO j=jj_begin,jj_end 413 DO i=ii_begin,ii_end 414 ij=(j-1)*iim+i 415 ps(ij)=psm1(ij)+tau*dps(ij) 416 ENDDO 404 !$SIMD 405 DO ij=ij_begin,ij_end 406 ps(ij)=psm1(ij)+tau*dps(ij) 417 407 ENDDO 418 408 ENDDO … … 428 418 IF (stage==1) THEN ! first stage : save model state in XXm1 429 419 DO l=ll_begin,ll_end 430 DO j=jj_begin,jj_end 431 DO i=ii_begin,ii_end 432 ij=(j-1)*iim+i 433 massm1(ij,l)=mass(ij,l) 434 ENDDO 435 ENDDO 420 !$SIMD 421 DO ij=ij_begin,ij_end 422 massm1(ij,l)=mass(ij,l) 423 ENDDO 436 424 ENDDO 437 425 END IF … … 439 427 ! updates are of the form x1 := x0 + tau*f(x1) 440 428 DO l=ll_begin,ll_end 441 DO j=jj_begin,jj_end 442 DO i=ii_begin,ii_end 443 ij=(j-1)*iim+i 444 mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l) 445 ENDDO 429 !$SIMD 430 DO ij=ij_begin,ij_end 431 mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l) 446 432 ENDDO 447 433 ENDDO … … 462 448 IF (stage==1) THEN ! first stage : save model state in XXm1 463 449 DO l=ll_begin,ll_end 464 DO j=jj_begin,jj_end 465 DO i=ii_begin,ii_end 466 ij=(j-1)*iim+i 450 !$SIMD 451 DO ij=ij_begin,ij_end 467 452 um1(ij+u_right,l)=u(ij+u_right,l) 468 453 um1(ij+u_lup,l)=u(ij+u_lup,l) 469 454 um1(ij+u_ldown,l)=u(ij+u_ldown,l) 470 455 theta_rhodzm1(ij,l)=theta_rhodz(ij,l) 471 ENDDO472 456 ENDDO 473 457 ENDDO … … 475 459 476 460 DO l=ll_begin,ll_end 477 DO j=jj_begin,jj_end 478 DO i=ii_begin,ii_end 479 ij=(j-1)*iim+i 461 !$SIMD 462 DO ij=ij_begin,ij_end 480 463 u(ij+u_right,l)=um1(ij+u_right,l)+tau*du(ij+u_right,l) 481 464 u(ij+u_lup,l)=um1(ij+u_lup,l)+tau*du(ij+u_lup,l) 482 465 u(ij+u_ldown,l)=um1(ij+u_ldown,l)+tau*du(ij+u_ldown,l) 483 466 theta_rhodz(ij,l)=theta_rhodzm1(ij,l)+tau*dtheta_rhodz(ij,l) 484 ENDDO485 467 ENDDO 486 468 ENDDO … … 564 546 565 547 DO l=ll_begin,ll_end 566 DO j=jj_begin-1,jj_end+1 567 DO i=ii_begin-1,ii_end+1 568 ij=(j-1)*iim+i 548 !$SIMD 549 DO ij=ij_begin_ext,ij_end_ext 569 550 hfluxt(ij+u_right,l) = tau*hflux(ij+u_right,l) 570 551 hfluxt(ij+u_lup,l) = tau*hflux(ij+u_lup,l) 571 552 hfluxt(ij+u_ldown,l) = tau*hflux(ij+u_ldown,l) 572 ENDDO573 553 ENDDO 574 554 ENDDO … … 576 556 IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag 577 557 DO l=ll_begin,ll_endp1 578 DO j=jj_begin,jj_end 579 DO i=ii_begin,ii_end 580 ij=(j-1)*iim+i 581 wfluxt(ij,l) = tau*wflux(ij,l) 582 ENDDO 558 !$SIMD 559 DO ij=ij_begin,ij_end 560 wfluxt(ij,l) = tau*wflux(ij,l) 583 561 ENDDO 584 562 ENDDO … … 588 566 589 567 DO l=ll_begin,ll_end 590 DO j=jj_begin-1,jj_end+1 591 DO i=ii_begin-1,ii_end+1 592 ij=(j-1)*iim+i 568 !$SIMD 569 DO ij=ij_begin_ext,ij_end_ext 593 570 hfluxt(ij+u_right,l) = hfluxt(ij+u_right,l)+tau*hflux(ij+u_right,l) 594 571 hfluxt(ij+u_lup,l) = hfluxt(ij+u_lup,l)+tau*hflux(ij+u_lup,l) 595 572 hfluxt(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)+tau*hflux(ij+u_ldown,l) 596 ENDDO597 573 ENDDO 598 574 ENDDO … … 600 576 IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag 601 577 DO l=ll_begin,ll_endp1 602 DO j=jj_begin,jj_end 603 DO i=ii_begin,ii_end 604 ij=(j-1)*iim+i 578 !$SIMD 579 DO ij=ij_begin,ij_end 605 580 wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l) 606 ENDDO607 581 ENDDO 608 582 ENDDO … … 613 587 END SUBROUTINE accumulate_fluxes 614 588 615 FUNCTION maxval_i(p)616 USE icosa617 IMPLICIT NONE618 REAL(rstd), DIMENSION(iim*jjm) :: p619 REAL(rstd) :: maxval_i620 INTEGER :: j, ij621 622 maxval_i=p((jj_begin-1)*iim+ii_begin)623 624 DO j=jj_begin-1,jj_end+1625 ij=(j-1)*iim626 maxval_i = MAX(maxval_i, MAXVAL(p(ij+ii_begin:ij+ii_end)))627 END DO628 END FUNCTION maxval_i629 630 FUNCTION maxval_ik(p)631 USE icosa632 IMPLICIT NONE633 REAL(rstd) :: p(iim*jjm, llm)634 REAL(rstd) :: maxval_ik(llm)635 INTEGER :: l,j, ij636 637 DO l=1,llm638 maxval_ik(l)=p((jj_begin-1)*iim+ii_begin,l)639 DO j=jj_begin-1,jj_end+1640 ij=(j-1)*iim641 maxval_ik(l) = MAX(maxval_ik(l), MAXVAL(p(ij+ii_begin:ij+ii_end,l)))642 END DO643 END DO644 END FUNCTION maxval_ik589 ! FUNCTION maxval_i(p) 590 ! USE icosa 591 ! IMPLICIT NONE 592 ! REAL(rstd), DIMENSION(iim*jjm) :: p 593 ! REAL(rstd) :: maxval_i 594 ! INTEGER :: j, ij 595 ! 596 ! maxval_i=p((jj_begin-1)*iim+ii_begin) 597 ! 598 ! DO j=jj_begin-1,jj_end+1 599 ! ij=(j-1)*iim 600 ! maxval_i = MAX(maxval_i, MAXVAL(p(ij+ii_begin:ij+ii_end))) 601 ! END DO 602 ! END FUNCTION maxval_i 603 604 ! FUNCTION maxval_ik(p) 605 ! USE icosa 606 ! IMPLICIT NONE 607 ! REAL(rstd) :: p(iim*jjm, llm) 608 ! REAL(rstd) :: maxval_ik(llm) 609 ! INTEGER :: l,j, ij 610 ! 611 ! DO l=1,llm 612 ! maxval_ik(l)=p((jj_begin-1)*iim+ii_begin,l) 613 ! DO j=jj_begin-1,jj_end+1 614 ! ij=(j-1)*iim 615 ! maxval_ik(l) = MAX(maxval_ik(l), MAXVAL(p(ij+ii_begin:ij+ii_end,l))) 616 ! END DO 617 ! END DO 618 ! END FUNCTION maxval_ik 645 619 646 620 END MODULE timeloop_gcm_mod
Note: See TracChangeset
for help on using the changeset viewer.