Changeset 148 for codes/icosagcm/trunk/src/advect.f90
- Timestamp:
- 03/18/13 15:44:08 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/advect.f90
r141 r148 11 11 !========================================================================== 12 12 13 SUBROUTINE init_advect(normal,tangent )13 SUBROUTINE init_advect(normal,tangent,one_over_sqrt_leng) 14 14 IMPLICIT NONE 15 15 REAL(rstd),INTENT(OUT) :: normal(3*iim*jjm,3) 16 16 REAL(rstd),INTENT(OUT) :: tangent(3*iim*jjm,3) 17 REAL(rstd),INTENT(OUT) :: one_over_sqrt_leng(iim*jjm) 17 18 18 19 INTEGER :: i,j,n … … 39 40 tangent(n+u_ldown,:)=xyz_v(n+z_down,:)-xyz_v(n+z_ldown,:) 40 41 tangent(n+u_ldown,:)=tangent(n+u_ldown,:)/sqrt(sum(tangent(n+u_ldown,:)**2)+1e-50) 41 END DO 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 42 47 ENDDO 43 48 … … 46 51 !======================================================================================= 47 52 48 SUBROUTINE compute_gradq3d(qi,gradq3d) 53 SUBROUTINE compute_gradq3d(qi,one_over_sqrt_leng,gradq3d) 54 USE trace 49 55 IMPLICIT NONE 50 56 REAL(rstd),INTENT(IN) :: qi(iim*jjm,llm) 57 REAL(rstd),INTENT(IN) :: one_over_sqrt_leng(iim*jjm) 51 58 REAL(rstd),INTENT(OUT) :: gradq3d(iim*jjm,llm,3) 52 59 REAL(rstd) :: maxq,minq,minq_c,maxq_c 53 REAL(rstd) :: alphamx,alphami,alpha ,maggrd ,leng60 REAL(rstd) :: alphamx,alphami,alpha ,maggrd 54 61 REAL(rstd) :: leng1,leng2 55 62 REAL(rstd) :: arr(2*iim*jjm) 63 REAL(rstd) :: ar(iim*jjm) 56 64 REAL(rstd) :: gradtri(2*iim*jjm,llm,3) 57 REAL(rstd) :: ar58 65 INTEGER :: i,j,n,k,ind,l 66 67 CALL trace_start("compute_gradq3d") 59 68 60 69 ! TODO : precompute ar, drop arr as output argument of gradq ? … … 67 76 DO i=ii_begin-1,ii_end+1 68 77 n=(j-1)*iim+i 69 CALL gradq(n, n+t_rup,n+t_lup,n+z_up,qi(:,l),gradtri(n+z_up,l,:),arr(n+z_up))70 CALL gradq(n, n+t_ldown,n+t_rdown,n+z_down,qi(:,l),gradtri(n+z_down,l,:),arr(n+z_down))78 CALL gradq(n,l,n+t_rup,n+t_lup,n+z_up,qi,gradtri(n+z_up,l,:),arr(n+z_up)) 79 CALL gradq(n,l,n+t_ldown,n+t_rdown,n+z_down,qi,gradtri(n+z_down,l,:),arr(n+z_down)) 71 80 END DO 72 81 END DO 73 82 END DO 74 83 75 ! Do l =1,llm 76 DO j=jj_begin,jj_end 77 DO i=ii_begin,ii_end 78 n=(j-1)*iim+i 79 gradq3d(n,:,:) = gradtri(n+z_up,:,:) + gradtri(n+z_down,:,:) + & 80 gradtri(n+z_rup,:,:) + gradtri(n+z_ldown,:,:) + & 81 gradtri(n+z_lup,:,:)+ gradtri(n+z_rdown,:,:) 82 ar = arr(n+z_up)+arr(n+z_lup)+arr(n+z_ldown)+arr(n+z_down)+arr(n+z_rdown)+arr(n+z_rup) 83 gradq3d(n,:,:) = gradq3d(n,:,:)/(ar+1.e-50) 84 END DO 85 END DO 86 ! END DO 84 DO j=jj_begin-1,jj_end+1 85 DO i=ii_begin-1,ii_end+1 86 n=(j-1)*iim+i 87 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 88 ENDDO 89 ENDDO 90 91 DO k=1,3 92 DO l =1,llm 93 DO j=jj_begin,jj_end 94 DO i=ii_begin,ii_end 95 n=(j-1)*iim+i 96 gradq3d(n,l,k) = ( gradtri(n+z_up,l,k) + gradtri(n+z_down,l,k) + & 97 gradtri(n+z_rup,l,k) + gradtri(n+z_ldown,l,k) + & 98 gradtri(n+z_lup,l,k)+ gradtri(n+z_rdown,l,k) ) / ar(n) 99 END DO 100 END DO 101 END DO 102 ENDDO 87 103 88 104 !============================================================================================= LIMITING … … 94 110 maggrd = dot_product(gradq3d(n,l,:),gradq3d(n,l,:)) 95 111 maggrd = sqrt(maggrd) 96 leng = max(sum((xyz_v(n+z_up,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_down,:) - xyz_i(n,:))**2), & 97 sum((xyz_v(n+z_rup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_rdown,:) - xyz_i(n,:))**2), & 98 sum((xyz_v(n+z_lup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_ldown,:) - xyz_i(n,:))**2)) 99 maxq_c = qi(n,l) + maggrd*sqrt(leng) 100 minq_c = qi(n,l) - maggrd*sqrt(leng) 112 ! leng = max(sum((xyz_v(n+z_up,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_down,:) - xyz_i(n,:))**2), & 113 ! sum((xyz_v(n+z_rup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_rdown,:) - xyz_i(n,:))**2), & 114 ! sum((xyz_v(n+z_lup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_ldown,:) - xyz_i(n,:))**2)) 115 ! maxq_c = qi(n,l) + maggrd*sqrt(leng(n)) 116 ! minq_c = qi(n,l) - maggrd*sqrt(leng(n)) 117 maxq_c = qi(n,l) + maggrd*one_over_sqrt_leng(n) 118 minq_c = qi(n,l) - maggrd*one_over_sqrt_leng(n) 101 119 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), & 102 120 qi(n+t_rdown,l),qi(n+t_ldown,l)) … … 112 130 END DO 113 131 END DO 132 133 CALL trace_end("compute_gradq3d") 114 134 END SUBROUTINE compute_gradq3d 115 135 116 136 ! Backward trajectories, for use with Miura approach 117 137 SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc) 138 USE trace 118 139 IMPLICIT NONE 119 140 REAL(rstd),INTENT(IN) :: normal(3*iim*jjm,3) … … 125 146 REAL(rstd) :: v_e(3), up_e, qe, ed(3) 126 147 INTEGER :: i,j,n,l 148 149 CALL trace_start("compute_backward_traj") 127 150 128 151 ! TODO : compute normal displacement ue*tau as hfluxt / mass(upwind) then reconstruct tangential displacement … … 181 204 ENDDO 182 205 END DO 206 207 CALL trace_end("compute_backward_traj") 208 183 209 END SUBROUTINE compute_backward_traj 184 210 … … 186 212 ! Slope-limited van Leer approach with hexagons 187 213 SUBROUTINE compute_advect_horiz(update_mass,hfluxt,cc,gradq3d, mass,qi) 214 USE trace 188 215 IMPLICIT NONE 189 216 LOGICAL, INTENT(IN) :: update_mass … … 197 224 REAL(rstd) :: qflux(3*iim*jjm,llm) 198 225 INTEGER :: i,j,n,k,l 226 227 CALL trace_start("compute_advect_horiz") 199 228 200 229 ! evaluate tracer field at point cc using piecewise linear reconstruction … … 262 291 END DO 263 292 293 CALL trace_end("compute_advect_horiz") 264 294 CONTAINS 265 295 … … 275 305 276 306 !========================================================================== 277 PURE SUBROUTINE gradq(n0, n1,n2,n3,q,dq,det)278 IMPLICIT NONE 279 INTEGER, INTENT(IN) :: n0, n1,n2,n3280 REAL(rstd), INTENT(IN) :: q(iim*jjm )307 PURE SUBROUTINE gradq(n0,l,n1,n2,n3,q,dq,det) 308 IMPLICIT NONE 309 INTEGER, INTENT(IN) :: n0,l,n1,n2,n3 310 REAL(rstd), INTENT(IN) :: q(iim*jjm,llm) 281 311 REAL(rstd), INTENT(OUT) :: dq(3), det 282 312 … … 293 323 A(3,1)=xyz_v(n3,1); A(3,2)= xyz_v(n3,2); A(3,3)=xyz_v(n3,3) 294 324 295 dq(1) = q(n1 )-q(n0)296 dq(2) = q(n2 )-q(n0)325 dq(1) = q(n1,l)-q(n0,l) 326 dq(2) = q(n2,l)-q(n0,l) 297 327 dq(3) = 0.0 298 328 ! CALL DGESV(3,1,A,3,IPIV,dq(:),3,info)
Note: See TracChangeset
for help on using the changeset viewer.