Changeset 138 for codes/icosagcm/trunk
- Timestamp:
- 02/16/13 17:03:57 (11 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/advect.f90
r137 r138 3 3 IMPLICIT NONE 4 4 5 PRIVATE 6 7 PUBLIC :: init_advect, compute_backward_traj, compute_gradq3d, compute_advect_horiz 8 5 9 CONTAINS 6 10 … … 8 12 9 13 SUBROUTINE init_advect(normal,tangent) 10 USE domain_mod11 USE dimensions12 USE geometry13 USE metric14 USE vector15 14 IMPLICIT NONE 16 15 REAL(rstd),INTENT(OUT) :: normal(3*iim*jjm,3) … … 48 47 49 48 SUBROUTINE compute_gradq3d(qi,gradq3d) 50 USE domain_mod51 USE dimensions52 USE geometry53 USE metric54 49 IMPLICIT NONE 55 50 REAL(rstd),INTENT(IN) :: qi(iim*jjm,llm) … … 62 57 REAL(rstd) :: ar 63 58 INTEGER :: i,j,n,k,ind,l 59 60 ! TODO : precompute ar, drop arr as output argument of gradq ? 61 64 62 !========================================================================================== GRADIENT 63 ! Compute gradient at triangles solving a linear system 64 ! arr = area of triangle joining centroids of hexagons 65 65 Do l = 1,llm 66 66 DO j=jj_begin-1,jj_end+1 … … 114 114 END SUBROUTINE compute_gradq3d 115 115 116 ! Backward trajectories, used inMiura approach116 ! Backward trajectories, for use with Miura approach 117 117 SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc) 118 USE domain_mod119 USE dimensions120 USE geometry121 USE metric122 118 IMPLICIT NONE 123 119 REAL(rstd),INTENT(IN) :: normal(3*iim*jjm,3) … … 127 123 REAL(rstd),INTENT(IN) :: tau 128 124 129 REAL(rstd) :: v_e(3), up_e 130 REAL(rstd) :: qe 131 REAL(rstd) :: ed(3) 125 REAL(rstd) :: v_e(3), up_e, qe, ed(3) 132 126 INTEGER :: i,j,n,l 133 127 128 ! TODO : compute normal displacement ue*tau as hfluxt / mass(upwind) then reconstruct tangential displacement 129 134 130 ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau 135 131 DO l = 1,llm … … 189 185 ! Horizontal transport (S. Dubey, T. Dubos) 190 186 ! Slope-limited van Leer approach with hexagons 191 SUBROUTINE compute_advect_horiz(update_mass, normal,tangent,qi,gradq3d,him,ue,hfluxt,bigt) 192 USE domain_mod 193 USE dimensions 194 USE geometry 195 USE metric 196 IMPLICIT NONE 197 LOGICAL, INTENT(IN) :: update_mass 198 REAL(rstd),INTENT(IN) :: normal(3*iim*jjm,3) 199 REAL(rstd),INTENT(IN) :: tangent(3*iim*jjm,3) 200 REAL(rstd),INTENT(INOUT) :: qi(iim*jjm,llm) 201 REAL(rstd),INTENT(IN) :: gradq3d(iim*jjm,llm,3) 202 REAL(rstd),INTENT(INOUT) :: him(iim*jjm,llm) 203 REAL(rstd),INTENT(IN) :: ue(iim*3*jjm,llm) 204 REAL(rstd),INTENT(IN) :: hfluxt(3*iim*jjm,llm) ! mass flux 205 REAL(rstd),INTENT(IN) :: bigt 206 207 REAL(rstd) :: dqi(iim*jjm,llm),dhi(iim*jjm,llm) 208 REAL(rstd) :: cc(3*iim*jjm,llm,3) ! barycenter of quadrilateral, where q is evaluated (1-point quadrature) 209 REAL(rstd) :: v_e(3*iim*jjm,llm,3) 210 REAL(rstd) :: up_e 211 212 ! REAL(rstd) :: qe(3*iim*jjm,llm) 213 REAL(rstd) :: qe 214 REAL(rstd) :: ed(3) 215 REAL(rstd) :: flxx(3*iim*jjm,llm) 187 SUBROUTINE compute_advect_horiz(update_mass,hfluxt,cc,gradq3d, mass,qi) 188 IMPLICIT NONE 189 LOGICAL, INTENT(IN) :: update_mass 190 REAL(rstd), INTENT(IN) :: gradq3d(iim*jjm,llm,3) 191 REAL(rstd), INTENT(IN) :: hfluxt(3*iim*jjm,llm) ! mass flux 192 REAL(rstd), INTENT(IN) :: cc(3*iim*jjm,llm,3) ! barycenter of quadrilateral, where q is evaluated (1-point quadrature) 193 REAL(rstd), INTENT(INOUT) :: mass(iim*jjm,llm) 194 REAL(rstd), INTENT(INOUT) :: qi(iim*jjm,llm) 195 196 REAL(rstd) :: dq,dmass,qe,ed(3), newmass 197 REAL(rstd) :: qflux(3*iim*jjm,llm) 216 198 INTEGER :: i,j,n,k,l 217 REAL(rstd):: him_old(llm) 218 219 220 ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*dt/2 221 ! TODO : this does not depend on q and should be done only once, not for each tracer 199 200 ! evaluate tracer field at point cc using piecewise linear reconstruction 201 ! q(cc)= q0 + gradq.(cc-xyz_i) with xi centroid of hexagon 202 ! ne*hfluxt>0 iff outgoing 222 203 DO l = 1,llm 223 204 DO j=jj_begin-1,jj_end+1 224 205 DO i=ii_begin-1,ii_end+1 225 206 n=(j-1)*iim+i 226 227 up_e =1/de(n+u_right)*( & 228 wee(n+u_right,1,1)*le(n+u_rup)*ue(n+u_rup,l)+ & 229 wee(n+u_right,2,1)*le(n+u_lup)*ue(n+u_lup,l)+ & 230 wee(n+u_right,3,1)*le(n+u_left)*ue(n+u_left,l)+ & 231 wee(n+u_right,4,1)*le(n+u_ldown)*ue(n+u_ldown,l)+ & 232 wee(n+u_right,5,1)*le(n+u_rdown)*ue(n+u_rdown,l)+ & 233 wee(n+u_right,1,2)*le(n+t_right+u_ldown)*ue(n+t_right+u_ldown,l)+ & 234 wee(n+u_right,2,2)*le(n+t_right+u_rdown)*ue(n+t_right+u_rdown,l)+ & 235 wee(n+u_right,3,2)*le(n+t_right+u_right)*ue(n+t_right+u_right,l)+ & 236 wee(n+u_right,4,2)*le(n+t_right+u_rup)*ue(n+t_right+u_rup,l)+ & 237 wee(n+u_right,5,2)*le(n+t_right+u_lup)*ue(n+t_right+u_lup,l) & 238 ) 239 v_e(n+u_right,l,:)= ue(n+u_right,l)*normal(n+u_right,:) + up_e*tangent(n+u_right,:) 240 241 up_e=1/de(n+u_lup)*( & 242 wee(n+u_lup,1,1)*le(n+u_left)*ue(n+u_left,l)+ & 243 wee(n+u_lup,2,1)*le(n+u_ldown)*ue(n+u_ldown,l)+ & 244 wee(n+u_lup,3,1)*le(n+u_rdown)*ue(n+u_rdown,l)+ & 245 wee(n+u_lup,4,1)*le(n+u_right)*ue(n+u_right,l)+ & 246 wee(n+u_lup,5,1)*le(n+u_rup)*ue(n+u_rup,l)+ & 247 wee(n+u_lup,1,2)*le(n+t_lup+u_right)*ue(n+t_lup+u_right,l)+ & 248 wee(n+u_lup,2,2)*le(n+t_lup+u_rup)*ue(n+t_lup+u_rup,l)+ & 249 wee(n+u_lup,3,2)*le(n+t_lup+u_lup)*ue(n+t_lup+u_lup,l)+ & 250 wee(n+u_lup,4,2)*le(n+t_lup+u_left)*ue(n+t_lup+u_left,l)+ & 251 wee(n+u_lup,5,2)*le(n+t_lup+u_ldown)*ue(n+t_lup+u_ldown,l) & 252 ) 253 v_e(n+u_lup,l,:)= ue(n+u_lup,l)*normal(n+u_lup,:) + up_e*tangent(n+u_lup,:) 254 255 up_e=1/de(n+u_ldown)*( & 256 wee(n+u_ldown,1,1)*le(n+u_rdown)*ue(n+u_rdown,l)+ & 257 wee(n+u_ldown,2,1)*le(n+u_right)*ue(n+u_right,l)+ & 258 wee(n+u_ldown,3,1)*le(n+u_rup)*ue(n+u_rup,l)+ & 259 wee(n+u_ldown,4,1)*le(n+u_lup)*ue(n+u_lup,l)+ & 260 wee(n+u_ldown,5,1)*le(n+u_left)*ue(n+u_left,l)+ & 261 wee(n+u_ldown,1,2)*le(n+t_ldown+u_lup)*ue(n+t_ldown+u_lup,l)+ & 262 wee(n+u_ldown,2,2)*le(n+t_ldown+u_left)*ue(n+t_ldown+u_left,l)+ & 263 wee(n+u_ldown,3,2)*le(n+t_ldown+u_ldown)*ue(n+t_ldown+u_ldown,l)+ & 264 wee(n+u_ldown,4,2)*le(n+t_ldown+u_rdown)*ue(n+t_ldown+u_rdown,l)+ & 265 wee(n+u_ldown,5,2)*le(n+t_ldown+u_right)*ue(n+t_ldown+u_right,l) & 266 ) 267 v_e(n+u_ldown,l,:)= ue(n+u_ldown,l)*normal(n+u_ldown,:) + up_e*tangent(n+u_ldown,:) 268 269 cc(n+u_right,l,:) = xyz_e(n+u_right,:) - v_e(n+u_right,l,:)*0.5*bigt !! redge is mid point of edge i 270 cc(n+u_lup,l,:) = xyz_e(n+u_lup,:) - v_e(n+u_lup,l,:)*0.5*bigt 271 cc(n+u_ldown,l,:) = xyz_e(n+u_ldown,:) - v_e(n+u_ldown,l,:)*0.5*bigt 272 ENDDO 273 ENDDO 274 END DO 275 ! end of tracer-independant computations 276 277 ! evaluate traver field at point cc usin piecewise linear reconstruction 278 ! q(cc)= q0 + gradq.(cc-xyz_i) with xi centroid of hexagon 279 DO l = 1,llm 280 DO j=jj_begin-1,jj_end+1 281 DO i=ii_begin-1,ii_end+1 282 n=(j-1)*iim+i 283 if (ne(n,right)*ue(n+u_right,l)>0) then 207 if (ne(n,right)*hfluxt(n+u_right,l)>0) then 284 208 ed = cc(n+u_right,l,:) - xyz_i(n,:) 285 qe =qi(n,l)+sum2(gradq3d(n,l,:),ed)209 qe = qi(n,l)+sum2(gradq3d(n,l,:),ed) 286 210 else 287 211 ed = cc(n+u_right,l,:) - xyz_i(n+t_right,:) 288 qe =qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed)212 qe = qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed) 289 213 endif 290 flxx(n+u_right,l) = hfluxt(n+u_right,:)*qe*ne(n,right)214 qflux(n+u_right,l) = hfluxt(n+u_right,l)*qe 291 215 292 if (ne(n,lup)* ue(n+u_lup,l)>0) then216 if (ne(n,lup)*hfluxt(n+u_lup,l)>0) then 293 217 ed = cc(n+u_lup,l,:) - xyz_i(n,:) 294 qe (n+u_lup,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed)218 qe = qi(n,l)+sum2(gradq3d(n,l,:),ed) 295 219 else 296 220 ed = cc(n+u_lup,l,:) - xyz_i(n+t_lup,:) 297 qe (n+u_lup,l)=qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed)221 qe = qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed) 298 222 endif 299 flxx(n+u_lup,:) = hfluxt(n+u_lup,:)*qe(n+u_lup,:)*ne(n,lup)223 qflux(n+u_lup,l) = hfluxt(n+u_lup,l)*qe 300 224 301 if (ne(n,ldown)* ue(n+u_ldown,l)>0) then225 if (ne(n,ldown)*hfluxt(n+u_ldown,l)>0) then 302 226 ed = cc(n+u_ldown,l,:) - xyz_i(n,:) 303 qe (n+u_ldown,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed)227 qe = qi(n,l)+sum2(gradq3d(n,l,:),ed) 304 228 else 305 229 ed = cc(n+u_ldown,l,:) - xyz_i(n+t_ldown,:) 306 qe (n+u_ldown,l)=qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed)230 qe = qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed) 307 231 endif 308 flxx(n+u_ldown,:) = hfluxt(n+u_ldown,:)*qe(n+u_ldown,:)*ne(n,ldown)232 qflux(n+u_ldown,l) = hfluxt(n+u_ldown,l)*qe 309 233 end do 310 234 end do 311 235 END DO 312 236 313 ! evaluate q flux as mass flux * qe 314 ! TODO : loop over k ! 315 ! TODO : merge with previous loop and eliminate unnecessary temporary qe 316 DO j=jj_begin-1,jj_end+1 317 DO i=ii_begin-1,ii_end+1 318 n=(j-1)*iim+i 319 ENDDO 320 ENDDO 321 322 ! update q and, if update_mass, mass 323 ! TODO : loop over k ! 324 DO j=jj_begin,jj_end 325 DO i=ii_begin,ii_end 326 n=(j-1)*iim+i 327 328 dhi(n,:)= -(1/Ai(n))*(he(n+u_right,:)*ne(n,right) + he(n+u_lup,:)*ne(n,lup) & 329 + he(n+u_ldown,:)*ne(n,ldown) + he(n+u_rup,:)*ne(n,rup) & 330 + he(n+u_left,:)*ne(n,left) + he(n+u_rdown,:)*ne(n,rdown) ) 331 332 dqi(n,:)= -(1/Ai(n))*(flxx(n+u_right,:)+flxx(n+u_lup,:) & 333 +flxx(n+u_ldown,:) - flxx(n+u_rup,:) & 334 - flxx(n+u_left,:) - flxx(n+u_rdown,:) ) 335 him_old(:) = him(n,:) 336 him(n,:) = him(n,:) + dhi(n,:) 337 qi(n,:) = (qi(n,:)*him_old(:) + dqi(n,:))/him(n,:) 338 237 ! update q and, if update_mass, update mass 238 DO l = 1,llm 239 DO j=jj_begin,jj_end 240 DO i=ii_begin,ii_end 241 n=(j-1)*iim+i 242 ! sign convention as Ringler et al. (2010) eq. 21 243 dmass = hfluxt(n+u_right,l)*ne(n,right) & 244 + hfluxt(n+u_lup,l) *ne(n,lup) & 245 + hfluxt(n+u_ldown,l)*ne(n,ldown) & 246 + hfluxt(n+u_rup,l) *ne(n,rup) & 247 + hfluxt(n+u_left,l) *ne(n,left) & 248 + hfluxt(n+u_rdown,l)*ne(n,rdown) 249 250 dq = qflux(n+u_right,l) *ne(n,right) & 251 + qflux(n+u_lup,l) *ne(n,lup) & 252 + qflux(n+u_ldown,l) *ne(n,ldown) & 253 + qflux(n+u_rup,l) *ne(n,rup) & 254 + qflux(n+u_left,l) *ne(n,left) & 255 + qflux(n+u_rdown,l) *ne(n,rdown) 256 257 newmass = mass(n,l) - dmass/Ai(n) 258 qi(n,l) = (qi(n,l)*mass(n,l) - dq/Ai(n) ) / newmass 259 IF(update_mass) mass(n,l) = newmass 260 END DO 339 261 END DO 340 262 END DO 341 263 342 264 CONTAINS 265 343 266 !==================================================================================== 344 REALFUNCTION sum2(a1,a2)267 PURE REAL(rstd) FUNCTION sum2(a1,a2) 345 268 IMPLICIT NONE 346 REAL ,INTENT(IN):: a1(3), a2(3)269 REAL(rstd),INTENT(IN):: a1(3), a2(3) 347 270 sum2 = a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3) 271 ! sum2 = 0. ! Godunov scheme 348 272 END FUNCTION sum2 349 273 350 274 END SUBROUTINE compute_advect_horiz 275 351 276 !========================================================================== 352 SUBROUTINE gradq(n0,n1,n2,n3,q,dq,det) 353 USE geometry 354 USE metric 355 USE dimensions 277 PURE SUBROUTINE gradq(n0,n1,n2,n3,q,dq,det) 356 278 IMPLICIT NONE 357 279 INTEGER, INTENT(IN) :: n0,n1,n2,n3 358 REAL,INTENT(IN) :: q(iim*jjm) 359 REAL,INTENT(OUT) :: dq(3) 360 REAL(rstd) ::det,detx,dety,detz 280 REAL(rstd), INTENT(IN) :: q(iim*jjm) 281 REAL(rstd), INTENT(OUT) :: dq(3), det 282 283 REAL(rstd) ::detx,dety,detz 361 284 INTEGER :: info 362 285 INTEGER :: IPIV(3) … … 365 288 REAL(rstd) :: B(3) 366 289 367 A(1,1)=xyz_i(n1,1) -xyz_i(n0,1); A(1,2)=xyz_i(n1,2)- xyz_i(n0,2); A(1,3)=xyz_i(n1,3) - xyz_i(n0,3)368 A( 2,1)=xyz_i(n2,1) - xyz_i(n0,1); A(2,2)=xyz_i(n2,2) - xyz_i(n0,2); A(2,3)=xyz_i(n2,3) -xyz_i(n0,3)369 A( 3,1)=xyz_v(n3,1); A(3,2)= xyz_v(n3,2); A(3,3)= xyz_v(n3,3)370 290 ! TODO : replace A by A1,A2,A3 291 A(1,1)=xyz_i(n1,1)-xyz_i(n0,1); A(1,2)=xyz_i(n1,2)-xyz_i(n0,2); A(1,3)=xyz_i(n1,3)-xyz_i(n0,3) 292 A(2,1)=xyz_i(n2,1)-xyz_i(n0,1); A(2,2)=xyz_i(n2,2)-xyz_i(n0,2); A(2,3)=xyz_i(n2,3)-xyz_i(n0,3) 293 A(3,1)=xyz_v(n3,1); A(3,2)= xyz_v(n3,2); A(3,3)=xyz_v(n3,3) 371 294 372 295 dq(1) = q(n1)-q(n0) … … 382 305 dq(3) = detz 383 306 END SUBROUTINE gradq 307 384 308 !========================================================================== 385 SUBROUTINE determinant(a1,a2,a3,det) 386 IMPLICIT NONE 387 REAL, DIMENSION(3) :: a1, a2,a3 388 REAL :: det,x1,x2,x3 309 PURE SUBROUTINE determinant(a1,a2,a3,det) 310 IMPLICIT NONE 311 REAL(rstd), DIMENSION(3), INTENT(IN) :: a1,a2,a3 312 REAL(rstd), INTENT(OUT) :: det 313 REAL(rstd) :: x1,x2,x3 389 314 x1 = a1(1) * (a2(2) * a3(3) - a2(3) * a3(2)) 390 315 x2 = a1(2) * (a2(1) * a3(3) - a2(3) * a3(1)) -
codes/icosagcm/trunk/src/advect_tracer.f90
r137 r138 1 1 MODULE advect_tracer_mod 2 2 USE icosa 3 IMPLICIT NONE 3 4 PRIVATE 4 5 … … 8 9 TYPE(t_field),POINTER :: f_cc(:) ! starting point of backward-trajectory (Miura approach) 9 10 10 ! TYPE(t_field),POINTER :: f_rhodzm1(:)11 ! TYPE(t_field),POINTER :: f_rhodz(:)12 13 11 REAL(rstd), PARAMETER :: pente_max=2.0 ! for vlz 14 12 … … 19 17 SUBROUTINE init_advect_tracer 20 18 USE advect_mod 21 IMPLICIT NONE22 19 REAL(rstd),POINTER :: tangent(:,:) 23 20 REAL(rstd),POINTER :: normal(:,:) 24 21 INTEGER :: ind 25 22 26 CALL allocate_field(f_normal,field_u,type_real,3) 27 CALL allocate_field(f_tangent,field_u,type_real,3) 28 CALL allocate_field(f_gradq3d,field_t,type_real,llm,3) 29 CALL allocate_field(f_cc,field_u,type_real,llm,3) 30 31 ! CALL allocate_field(f_rhodz,field_t,type_real,llm) 32 ! CALL allocate_field(f_rhodzm1,field_t,type_real,llm) 23 CALL allocate_field(f_normal,field_u,type_real,3, name='normal') 24 CALL allocate_field(f_tangent,field_u,type_real,3, name='tangent') 25 CALL allocate_field(f_gradq3d,field_t,type_real,llm,3, name='gradq3d') 26 CALL allocate_field(f_cc,field_u,type_real,llm,3, name='cc') 33 27 34 28 DO ind=1,ndomain … … 43 37 44 38 SUBROUTINE advect_tracer(f_hfluxt, f_wfluxt,f_u, f_q,f_rhodz) 45 USE icosa46 USE field_mod47 USE domain_mod48 USE dimensions49 USE grid_param50 USE geometry51 USE metric52 39 USE advect_mod 53 USE advect_mod54 USE disvert_mod55 40 USE mpipara 56 41 IMPLICIT NONE … … 62 47 TYPE(t_field),POINTER :: f_rhodz(:) ! mass field at beginning of macro time step 63 48 64 REAL(rstd),POINTER :: q(:,:,:), normal(:,: ,:), tangent(:,:,:), gradq3d(:,:,:)49 REAL(rstd),POINTER :: q(:,:,:), normal(:,:), tangent(:,:), gradq3d(:,:,:), cc(:,:,:) 65 50 REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:) 66 51 REAL(rstd),POINTER :: rhodz(:,:), u(:,:) 67 INTEGER :: ind 52 INTEGER :: ind,k 68 53 69 54 CALL transfert_request(f_u,req_e1) 70 CALL transfert_request(f_u,req_e1) 55 ! CALL transfert_request(f_hfluxt,req_e1) ! BUG : This (unnecessary) transfer makes the computation go wrong 56 CALL transfert_request(f_wfluxt,req_i1) 71 57 CALL transfert_request(f_q,req_i1) 58 CALL transfert_request(f_rhodz,req_i1) 72 59 73 60 IF (is_mpi_root) PRINT *, 'Advection scheme' 74 61 62 ! DO ind=1,ndomain 63 ! CALL swap_dimensions(ind) 64 ! CALL swap_geometry(ind) 65 ! normal = f_normal(ind) 66 ! tangent = f_tangent(ind) 67 ! cc = f_cc(ind) 68 ! u = f_u(ind) 69 ! q = f_q(ind) 70 ! rhodz = f_rhodz(ind) 71 ! hfluxt = f_hfluxt(ind) 72 ! wfluxt = f_wfluxt(ind) 73 ! gradq3d = f_gradq3d(ind) 74 ! 75 ! ! 1/2 vertical transport 76 ! DO k = 1, nqtot 77 ! CALL vlz(k==nqtot,0.5, wfluxt,rhodz,q(:,:,k)) 78 ! END DO 79 ! 80 ! ! horizontal transport 81 ! CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc) 82 ! DO k = 1,nqtot 83 ! CALL compute_gradq3d(q(:,:,k),gradq3d) 84 ! CALL compute_advect_horiz(k==nqtot,hfluxt,cc,gradq3d, rhodz,q(:,:,k)) 85 ! END DO 86 ! 87 ! ! 1/2 vertical transport 88 ! DO k = 1,nqtot 89 ! CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k)) 90 ! END DO 91 ! END DO 92 93 ! 1/2 vertical transport + back-trajectories 94 DO ind=1,ndomain 95 CALL swap_dimensions(ind) 96 CALL swap_geometry(ind) 97 normal = f_normal(ind) 98 tangent = f_tangent(ind) 99 cc = f_cc(ind) 100 u = f_u(ind) 101 q = f_q(ind) 102 rhodz = f_rhodz(ind) 103 wfluxt = f_wfluxt(ind) 104 DO k = 1, nqtot 105 CALL vlz(k==nqtot,0.5, wfluxt,rhodz,q(:,:,k)) 106 END DO 107 CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc) 108 END DO 109 110 CALL transfert_request(f_q,req_i1) ! necessary ? 111 CALL transfert_request(f_rhodz,req_i1) ! necessary ? 112 113 ! horizontal transport - split in two to place transfer of gradq3d 114 115 DO k = 1, nqtot 116 DO ind=1,ndomain 117 CALL swap_dimensions(ind) 118 CALL swap_geometry(ind) 119 q = f_q(ind) 120 gradq3d = f_gradq3d(ind) 121 CALL compute_gradq3d(q(:,:,k),gradq3d) 122 END DO 123 124 CALL transfert_request(f_gradq3d,req_i1) 125 126 DO ind=1,ndomain 127 CALL swap_dimensions(ind) 128 CALL swap_geometry(ind) 129 cc = f_cc(ind) 130 q = f_q(ind) 131 rhodz = f_rhodz(ind) 132 hfluxt = f_hfluxt(ind) 133 gradq3d = f_gradq3d(ind) 134 CALL compute_advect_horiz(k==nqtot,hfluxt,cc,gradq3d, rhodz,q(:,:,k)) 135 END DO 136 END DO 137 138 CALL transfert_request(f_q,req_i1) ! necessary ? 139 CALL transfert_request(f_rhodz,req_i1) ! necessary ? 140 141 ! 1/2 vertical transport 75 142 DO ind=1,ndomain 76 143 CALL swap_dimensions(ind) … … 78 145 q = f_q(ind) 79 146 rhodz = f_rhodz(ind) 80 hfluxt = f_hfluxt(ind)81 147 wfluxt = f_wfluxt(ind) 82 normal = f_normal(ind) 83 tangent = f_tangent(ind) 84 gradq3d = f_gradq3d(ind) 85 CALL compute_adv_tracer(tangent,normal,hfluxt,wfluxt,u, q,rhodz,gradq3d) 148 DO k = 1,nqtot 149 CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k)) 150 END DO 86 151 END DO 87 152 88 153 END SUBROUTINE advect_tracer 89 154 90 SUBROUTINE compute_adv_tracer(tangent,normal,hfluxt,wfluxt,u, q,rhodz,gradq3d)91 USE icosa92 USE field_mod93 USE domain_mod94 USE dimensions95 USE grid_param96 USE geometry97 USE metric98 USE advect_mod99 USE advect_mod100 USE disvert_mod101 USE mpipara102 REAL(rstd), INTENT(IN) :: tangent(:,:,:), normal(:,:,:)103 REAL(rstd), INTENT(IN) :: hfluxt(:,:), wfluxt(:,:), u(:,:)104 REAL(rstd), INTENT(INOUT) :: rhodz(:,:)105 REAL(rstd), INTENT(INOUT) :: q(:,:,:)106 REAL(rstd), INTENT(OUT) :: gradq3d(:,:,:)107 INTEGER :: k108 ! for mass/tracer consistency each transport step also updates the mass field rhodz109 ! mass is updated after all tracers have been updated, when k==nqtot110 111 ! 1/2 vertical transport112 DO k = 1, nqtot113 CALL vlz(k==nqtot, 0.5,wfluxt,rhodz,q(:,:,k))114 END DO115 116 ! horizontal transport117 DO k = 1,nqtot118 CALL compute_gradq3d(q(:,:,k),gradq3d)119 CALL compute_advect_horiz(k==nqtot, tangent,normal,q(:,:,k),gradq3d,rhodz,u,hfluxt,dt*itau_adv)120 END DO121 122 ! 1/2 vertical transport123 DO k = 1,nqtot124 CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k))125 END DO126 END SUBROUTINE compute_adv_tracer127 128 155 SUBROUTINE vlz(update_mass, fac,wfluxt,mass, q) 129 156 ! … … 135 162 ! wfluxt >0 for upward transport 136 163 ! ******************************************************************** 137 USE icosa138 164 IMPLICIT NONE 139 165 LOGICAL, INTENT(IN) :: update_mass … … 152 178 ! finite difference of q 153 179 DO l=2,llm 154 DO j=jj_begin ,jj_end155 DO i=ii_begin ,ii_end180 DO j=jj_begin-1,jj_end+1 181 DO i=ii_begin-1,ii_end+1 156 182 ij=(j-1)*iim+i 157 183 dzqw(ij,l)=q(ij,l)-q(ij,l-1) … … 164 190 ! dzq = slope*dz, i.e. the reconstructed q varies by dzq inside level l 165 191 DO l=2,llm-1 166 DO j=jj_begin ,jj_end167 DO i=ii_begin ,ii_end192 DO j=jj_begin-1,jj_end+1 193 DO i=ii_begin-1,ii_end+1 168 194 ij=(j-1)*iim+i 169 195 IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN … … 179 205 180 206 ! 0 slope in top and bottom layers 181 DO j=jj_begin ,jj_end182 DO i=ii_begin ,ii_end207 DO j=jj_begin-1,jj_end+1 208 DO i=ii_begin-1,ii_end+1 183 209 ij=(j-1)*iim+i 184 210 dzq(ij,1)=0. … … 190 216 ! then amount of q leaving level l/l+1 = wq = w * qq 191 217 DO l = 1,llm-1 192 DO j=jj_begin,jj_end 193 DO i=ii_begin,ii_end 194 ij=(j-1)*iim+i 195 IF(wfluxt(ij,l+1).gt.0.) THEN 196 ! downward transport, sigw<0 197 ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0 198 w = fac*wfluxt(ij,l+1) 218 DO j=jj_begin-1,jj_end+1 219 DO i=ii_begin-1,ii_end+1 220 ij=(j-1)*iim+i 221 w = fac*wfluxt(ij,l+1) 222 IF(w>0.) THEN ! upward transport, upwind side is at level l 223 sigw = w/mass(ij,l) 224 qq = q(ij,l)+0.5*(1.-sigw)*dzq(ij,l) ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0 225 ELSE ! downward transport, upwind side is at level l+1 199 226 sigw = w/mass(ij,l+1) 200 qq = q(ij,l+1) - 0.5*(1.+sigw)*dzq(ij,l+1) 201 wq(ij,l+1) = w*qq 202 ELSE 203 ! upward transport, sigw>0 204 ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0 205 w = fac*wfluxt(ij,l) 206 sigw = w/mass(ij,l) 207 qq = q(ij,l)+0.5*(1.-sigw)*dzq(ij,l) 208 wq(ij,l+1) = w*qq 227 qq = q(ij,l+1)-0.5*(1.+sigw)*dzq(ij,l+1) ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0 209 228 ENDIF 229 wq(ij,l+1) = w*qq 210 230 ENDDO 211 231 ENDDO 212 232 END DO 213 233 ! wq = 0 at top and bottom 214 DO j=jj_begin ,jj_end215 DO i=ii_begin ,ii_end234 DO j=jj_begin-1,jj_end+1 235 DO i=ii_begin-1,ii_end+1 216 236 ij=(j-1)*iim+i 217 237 wq(ij,llm+1)=0. … … 222 242 ! update q, mass is updated only after all q's have been updated 223 243 DO l=1,llm 224 DO j=jj_begin ,jj_end225 DO i=ii_begin ,ii_end244 DO j=jj_begin-1,jj_end+1 245 DO i=ii_begin-1,ii_end+1 226 246 ij=(j-1)*iim+i 227 247 newmass = mass(ij,l) + fac*(wfluxt(ij,l)-wfluxt(ij,l+1)) -
codes/icosagcm/trunk/src/caldyn_gcm.f90
r134 r138 297 297 REAL(rstd),INTENT(IN) :: qu(iim*3*jjm,llm) 298 298 299 REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm), hflux(iim*3*jjm,llm) 299 REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm), hflux(iim*3*jjm,llm) ! hflux in kg/s 300 300 REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm) 301 301 REAL(rstd),INTENT(OUT) :: dps(iim*jjm) 302 REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux 302 REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s) 303 303 304 304 INTEGER :: i,j,ij,l … … 312 312 REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:) ! geopotential 313 313 REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux 314 REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:) ! mass flux convergence314 REAL(rstd),ALLOCATABLE,SAVE :: divm(:,:) ! mass flux divergence 315 315 REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:) ! Bernouilli function 316 316 … … 328 328 ALLOCATE(phi(iim*jjm,llm)) ! geopotential 329 329 ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux 330 ALLOCATE( convm(iim*jjm,llm)) ! mass flux convergence330 ALLOCATE(divm(iim*jjm,llm)) ! mass flux convergence 331 331 ALLOCATE(berni(iim*jjm,llm)) ! bernouilli term 332 332 … … 363 363 DO i=ii_begin,ii_end 364 364 ij=(j-1)*iim+i 365 ! convm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21366 convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*hflux(ij+u_right,l) + &365 ! divm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 366 divm(ij,l)= 1./Ai(ij)*(ne(ij,right)*hflux(ij+u_right,l) + & 367 367 ne(ij,rup)*hflux(ij+u_rup,l) + & 368 368 ne(ij,lup)*hflux(ij+u_lup,l) + & … … 383 383 ENDDO 384 384 385 !!! cumulate mass flux convergence from top to bottom385 !!! cumulate mass flux divergence from top to bottom 386 386 DO l = llm-1, 1, -1 387 387 !$OMP DO … … 389 389 DO i=ii_begin,ii_end 390 390 ij=(j-1)*iim+i 391 convm(ij,l) = convm(ij,l) + convm(ij,l+1)391 divm(ij,l) = divm(ij,l) + divm(ij,l+1) 392 392 ENDDO 393 393 ENDDO … … 402 402 ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt 403 403 ! => w>0 for upward transport 404 wflux( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )404 wflux( ij, l+1 ) = divm( ij, l+1 ) - bp(l+1) * divm( ij, 1 ) 405 405 ENDDO 406 406 ENDDO 407 407 ENDDO 408 408 409 ! compute dps, vertical mass flux at the surface =0409 ! compute dps, set vertical mass flux at the surface to 0 410 410 !$OMP DO 411 411 DO j=jj_begin,jj_end … … 414 414 wflux(ij,1) = 0. 415 415 ! dps/dt = -int(div flux)dz 416 dps(ij)=- convm(ij,1) * g416 dps(ij)=-divm(ij,1) * g 417 417 ENDDO 418 418 ENDDO … … 644 644 DEALLOCATE(phi) ! geopotential 645 645 DEALLOCATE(Ftheta) ! theta flux 646 DEALLOCATE( convm) ! mass flux convergence646 DEALLOCATE(divm) ! mass flux divergence 647 647 DEALLOCATE(berni) ! bernouilli term 648 648 !!$OMP END MASTER -
codes/icosagcm/trunk/src/etat0_dcmip4.f90
r80 r138 174 174 175 175 q(:,:,1)=theta(:,:) 176 IF(nqtot>2) q(:,:,3)=1. 176 177 177 178 DO l=1,llm … … 200 201 -cos(latc)*sin(lat)*cos(lon-lonc))/(sqrt(1-K**2))) 201 202 q(ij,l,2)=ABS(g/preff*(-1./radius*duodeta*dthetaodlat-(2*sin(lat)*omega+vort)*dthetaodeta)) 203 IF(nqtot>3) q(ij,l,4)=cos(lon)*cos(lat) 202 204 ENDDO 203 205 ENDDO -
codes/icosagcm/trunk/src/field.f90
r29 r138 11 11 12 12 TYPE t_field 13 CHARACTER(30) :: name 13 14 REAL(rstd),POINTER :: rval2d(:) 14 15 REAL(rstd),POINTER :: rval3d(:,:) … … 45 46 CONTAINS 46 47 47 SUBROUTINE allocate_field(field,field_type,data_type,dim1,dim2 )48 SUBROUTINE allocate_field(field,field_type,data_type,dim1,dim2,name) 48 49 USE domain_mod 49 50 IMPLICIT NONE … … 52 53 INTEGER,INTENT(IN) :: data_type 53 54 INTEGER,OPTIONAL :: dim1,dim2 55 CHARACTER(*), OPTIONAL :: name 54 56 INTEGER :: ind 55 57 INTEGER :: ii_size,jj_size 56 58 57 ALLOCATE(field(ndomain)) 59 ALLOCATE(field(ndomain)) 58 60 59 61 DO ind=1,ndomain 60 62 IF(PRESENT(name)) THEN 63 field(ind)%name = name 64 ELSE 65 field(ind)%name = '(unkown)' 66 END IF 67 61 68 IF (PRESENT(dim2)) THEN 62 69 field(ind)%ndim=4 … … 224 231 TYPE(t_field),INTENT(IN) :: field 225 232 226 IF (field%ndim/=2 .OR. field%data_type/=type_real) STOP 'get_val_r2d : bad pointer assignation' 233 IF (field%ndim/=2 .OR. field%data_type/=type_real) THEN 234 PRINT *, 'get_val_r2d : bad pointer assignation with ' // TRIM(field%name) 235 STOP 236 END IF 227 237 field_pt=>field%rval2d 228 238 END SUBROUTINE getval_r2d … … 233 243 TYPE(t_field),INTENT(IN) :: field 234 244 235 IF (field%ndim/=3 .OR. field%data_type/=type_real) STOP 'get_val_r3d : bad pointer assignation' 245 IF (field%ndim/=3 .OR. field%data_type/=type_real) THEN 246 PRINT *, 'get_val_r3d : bad pointer assignation with ' // TRIM(field%name) 247 STOP 248 END IF 236 249 field_pt=>field%rval3d 237 250 END SUBROUTINE getval_r3d … … 242 255 TYPE(t_field),INTENT(IN) :: field 243 256 244 IF (field%ndim/=4 .OR. field%data_type/=type_real) STOP 'get_val_r4d : bad pointer assignation' 257 IF (field%ndim/=4 .OR. field%data_type/=type_real) THEN 258 PRINT *, 'get_val_r4d : bad pointer assignation with ' // TRIM(field%name) 259 STOP 260 END IF 245 261 field_pt=>field%rval4d 246 END SUBROUTINE getval_r4d 247 262 END SUBROUTINE getval_r4d 248 263 249 264 -
codes/icosagcm/trunk/src/timeloop_gcm.f90
r135 r138 46 46 INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period, scheme 47 47 CHARACTER(len=255) :: scheme_name 48 LOGICAL :: fluxt_zero ! set to .TRUE. to start accumulating fluxes in time48 LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time 49 49 ! INTEGER :: itaumax 50 50 ! REAL(rstd) ::write_period … … 141 141 CALL swap_geometry(ind) 142 142 rhodz=f_rhodz(ind); ps=f_ps(ind) 143 CALL compute_rhodz( ps,rhodz) ! save rhodz for transport scheme before dynamics update ps143 CALL compute_rhodz(.TRUE., ps,rhodz) ! save rhodz for transport scheme before dynamics update ps 144 144 END DO 145 fluxt_zero=.FALSE. 146 145 fluxt_zero=.TRUE. 146 147 ! check that rhodz is consistent with ps 148 CALL transfert_request(f_rhodz,req_i1) 149 CALL transfert_request(f_ps,req_i1) 150 DO ind=1,ndomain 151 CALL swap_dimensions(ind) 152 CALL swap_geometry(ind) 153 rhodz=f_rhodz(ind); ps=f_ps(ind) 154 CALL compute_rhodz(.FALSE., ps, rhodz) 155 END DO 156 147 157 DO it=0,itaumax 148 158 … … 182 192 183 193 IF(MOD(it+1,itau_adv)==0) THEN 194 CALL transfert_request(f_wfluxt,req_i1) ! FIXME 195 ! CALL transfert_request(f_hfluxt,req_e1) ! FIXME 196 184 197 CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz) ! update q and rhodz after RK step 185 198 fluxt_zero=.TRUE. 199 200 ! FIXME : check that rhodz is consistent with ps 201 CALL transfert_request(f_rhodz,req_i1) 202 CALL transfert_request(f_ps,req_i1) 203 CALL transfert_request(f_dps,req_i1) ! FIXME 204 CALL transfert_request(f_wflux,req_i1) ! FIXME 205 DO ind=1,ndomain 206 CALL swap_dimensions(ind) 207 CALL swap_geometry(ind) 208 rhodz=f_rhodz(ind); ps=f_ps(ind); dps=f_dps(ind); 209 wflux=f_wflux(ind); wfluxt=f_wfluxt(ind) 210 CALL compute_rhodz(.FALSE., ps, rhodz) 211 END DO 212 186 213 END IF 187 214 … … 195 222 LOGICAL :: with_dps 196 223 INTEGER :: ind 197 198 224 DO ind=1,ndomain 225 CALL swap_dimensions(ind) 226 CALL swap_geometry(ind) 199 227 IF(with_dps) THEN 200 228 ps=f_ps(ind) ; dps=f_dps(ind) ; 201 229 ps(:)=ps(:)+dt*dps(:) 202 230 hflux=f_hflux(ind); hfluxt=f_hfluxt(ind) 203 wflux=f_ hflux(ind); wfluxt=f_wfluxt(ind)204 CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero )231 wflux=f_wflux(ind); wfluxt=f_wfluxt(ind) 232 CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind)) 205 233 END IF 206 234 u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) … … 221 249 222 250 DO ind=1,ndomain 251 CALL swap_dimensions(ind) 252 CALL swap_geometry(ind) 223 253 ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) 224 254 psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind) … … 234 264 IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage 235 265 hflux=f_hflux(ind); hfluxt=f_hfluxt(ind) 236 wflux=f_ hflux(ind); wfluxt=f_wfluxt(ind)237 CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt,dt,fluxt_zero)266 wflux=f_wflux(ind); wfluxt=f_wfluxt(ind) 267 CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind)) 238 268 END IF 239 269 END DO … … 245 275 246 276 DO ind=1,ndomain 277 CALL swap_dimensions(ind) 278 CALL swap_geometry(ind) 247 279 ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) 248 280 psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind) … … 277 309 tau = dt/nb_stage 278 310 DO ind=1,ndomain 311 CALL swap_dimensions(ind) 312 CALL swap_geometry(ind) 313 279 314 ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) 280 315 psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind) … … 320 355 321 356 DO ind=1,ndomain 357 CALL swap_dimensions(ind) 358 CALL swap_geometry(ind) 322 359 ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) 323 360 dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind) … … 343 380 END SUBROUTINE timeloop 344 381 345 SUBROUTINE compute_rhodz( ps, rhodz)382 SUBROUTINE compute_rhodz(comp, ps, rhodz) 346 383 USE icosa 347 384 USE disvert_mod 385 LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check 348 386 REAL(rstd), INTENT(IN) :: ps(iim*jjm) 349 387 REAL(rstd), INTENT(OUT) :: rhodz(iim*jjm,llm) 350 INTEGER :: l,i,j,ij 388 REAL(rstd) :: m, err 389 INTEGER :: l,i,j,ij,dd 390 err=0. 391 IF(comp) THEN 392 dd=1 393 ELSE 394 ! dd=-1 395 dd=0 396 END IF 397 351 398 DO l = 1, llm 352 DO j=jj_begin- 1,jj_end+1353 DO i=ii_begin- 1,ii_end+1399 DO j=jj_begin-dd,jj_end+dd 400 DO i=ii_begin-dd,ii_end+dd 354 401 ij=(j-1)*iim+i 355 rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g 402 m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 403 IF(comp) THEN 404 rhodz(ij,l) = m 405 ELSE 406 err = MAX(err,abs(m-rhodz(ij,l))) 407 END IF 356 408 ENDDO 357 409 ENDDO 358 410 ENDDO 411 412 IF(.NOT. comp) THEN 413 IF(err>1e-10) THEN 414 PRINT *, 'Discrepancy between ps and rhodz detected', err 415 STOP 416 ELSE 417 ! PRINT *, 'No discrepancy between ps and rhodz detected' 418 END IF 419 END IF 420 359 421 END SUBROUTINE compute_rhodz 360 422 361 SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt,tau,fluxt_zero)423 SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero) 362 424 USE icosa 363 425 REAL(rstd), INTENT(IN) :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1) … … 366 428 LOGICAL, INTENT(INOUT) :: fluxt_zero 367 429 IF(fluxt_zero) THEN 430 ! PRINT *, 'Accumulating fluxes (first)' 368 431 fluxt_zero=.FALSE. 369 432 hfluxt = tau*hflux 370 433 wfluxt = tau*wflux 371 434 ELSE 435 ! PRINT *, 'Accumulating fluxes (next)' 372 436 hfluxt = hfluxt + tau*hflux 373 437 wfluxt = wfluxt + tau*wflux
Note: See TracChangeset
for help on using the changeset viewer.