Changeset 136 for codes/icosagcm/trunk
- Timestamp:
- 02/15/13 13:42:05 (12 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/advect.f90
r69 r136 115 115 116 116 !=================================================================================================== 117 SUBROUTINE compute_advect_horiz( normal,tangent,qi,gradq3d,him,ue,he,bigt)117 SUBROUTINE compute_advect_horiz(update_mass, normal,tangent,qi,gradq3d,him,ue,he,bigt) 118 118 USE domain_mod 119 119 USE dimensions … … 121 121 USE metric 122 122 IMPLICIT NONE 123 LOGICAL, INTENT(IN) :: update_mass 123 124 REAL(rstd),INTENT(IN) :: normal(3*iim*jjm,3) 124 125 REAL(rstd),INTENT(IN) :: tangent(3*iim*jjm,3) … … 131 132 132 133 REAL(rstd) :: dqi(iim*jjm,llm),dhi(iim*jjm,llm) 133 REAL(rstd) :: cc(3*iim*jjm,llm,3) 134 REAL(rstd) :: cc(3*iim*jjm,llm,3) ! barycenter of quadrilateral, where q is evaluated (1-point quadrature) 134 135 REAL(rstd) :: v_e(3*iim*jjm,llm,3) 135 136 REAL(rstd) :: up_e … … 142 143 143 144 144 !go to 123 145 ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*dt/2 146 ! TODO : this does not depend on q and should be done only once, not for each tracer 145 147 DO l = 1,llm 146 148 DO j=jj_begin-1,jj_end+1 … … 196 198 ENDDO 197 199 END DO 198 !123 continue 199 !========================================================================================================== 200 ! end of tracer-independant computations 201 202 ! evaluate traver field at point cc usin piecewise linear reconstruction 203 ! q(cc)= q0 + gradq.(cc-xyz_i) with xi centroid of hexagon 200 204 DO l = 1,llm 201 205 DO j=jj_begin-1,jj_end+1 … … 227 231 END DO 228 232 229 233 ! evaluate q flux as mass flux * qe 234 ! TODO : loop over k ! 235 ! TODO : merge with previous loop and eliminate unnecessary temporary qe 230 236 DO j=jj_begin-1,jj_end+1 231 237 DO i=ii_begin-1,ii_end+1 … … 237 243 ENDDO 238 244 245 ! update q and, if update_mass, mass 246 ! TODO : loop over k ! 239 247 DO j=jj_begin,jj_end 240 248 DO i=ii_begin,ii_end -
codes/icosagcm/trunk/src/advect_tracer.f90
r133 r136 2 2 USE icosa 3 3 PRIVATE 4 INTEGER,PARAMETER::iapp_tracvl= 15 4 6 5 TYPE(t_field),POINTER :: f_normal(:) … … 8 7 TYPE(t_field),POINTER :: f_gradq3d(:) 9 8 10 PUBLIC init_advect_tracer, advect_tracer_rhodz, advect_tracer 9 TYPE(t_field),POINTER :: f_wg(:) 10 TYPE(t_field),POINTER :: f_zm(:) 11 TYPE(t_field),POINTER :: f_zq(:) 12 13 TYPE(t_field),POINTER :: f_massflxc(:) 14 TYPE(t_field),POINTER :: f_massflx(:) 15 TYPE(t_field),POINTER :: f_uc(:) 16 TYPE(t_field),POINTER :: f_rhodzm1(:) 17 TYPE(t_field),POINTER :: f_rhodz(:) 18 19 REAL(rstd), PARAMETER :: pente_max=2.0 ! for vlz 20 21 PUBLIC init_advect_tracer, advect_tracer 11 22 12 23 CONTAINS … … 22 33 CALL allocate_field(f_tangent,field_u,type_real,3) 23 34 CALL allocate_field(f_gradq3d,field_t,type_real,llm,3) 35 36 CALL allocate_field(f_wg,field_t,type_real,llm) 37 CALL allocate_field(f_zm,field_t,type_real,llm) 38 CALL allocate_field(f_zq,field_t,type_real,llm,nqtot) 39 40 CALL allocate_field(f_rhodz,field_t,type_real,llm) 41 CALL allocate_field(f_rhodzm1,field_t,type_real,llm) 42 CALL allocate_field(f_massflxc,field_u,type_real,llm) 43 CALL allocate_field(f_massflx,field_u,type_real,llm) 44 CALL allocate_field(f_uc,field_u,type_real,llm) 24 45 25 46 DO ind=1,ndomain … … 33 54 END SUBROUTINE init_advect_tracer 34 55 35 SUBROUTINE advect_tracer(f_ ps,f_u,f_q)56 SUBROUTINE advect_tracer(f_hfluxt, f_wfluxt,f_u, f_q,f_rhodz) 36 57 USE icosa 37 USE advect_mod38 USE disvert_mod39 USE mpipara40 IMPLICIT NONE41 TYPE(t_field),POINTER :: f_ps(:)42 TYPE(t_field),POINTER :: f_u(:)43 TYPE(t_field),POINTER :: f_q(:)44 REAL(rstd),POINTER :: q(:,:,:)45 REAL(rstd),POINTER :: u(:,:)46 REAL(rstd),POINTER :: ps(:)47 REAL(rstd),POINTER :: massflx(:,:)48 REAL(rstd),POINTER :: rhodz(:,:)49 TYPE(t_field),POINTER,SAVE :: f_massflxc(:)50 TYPE(t_field),POINTER,SAVE :: f_massflx(:)51 TYPE(t_field),POINTER,SAVE :: f_uc(:)52 TYPE(t_field),POINTER,SAVE :: f_rhodzm1(:)53 TYPE(t_field),POINTER,SAVE :: f_rhodz(:)54 REAL(rstd),POINTER,SAVE :: massflxc(:,:)55 REAL(rstd),POINTER,SAVE :: uc(:,:)56 REAL(rstd),POINTER,SAVE :: rhodzm1(:,:)57 REAL(rstd):: bigt58 INTEGER :: ind,it,i,j,l,ij59 INTEGER,SAVE :: iadvtr=060 LOGICAL,SAVE:: first=.TRUE.61 !------------------------------------------------------sarvesh62 CALL transfert_request(f_ps,req_i1)63 CALL transfert_request(f_u,req_e1)64 CALL transfert_request(f_u,req_e1)65 CALL transfert_request(f_q,req_i1)66 67 IF ( first ) THEN68 CALL allocate_field(f_rhodz,field_t,type_real,llm)69 CALL allocate_field(f_rhodzm1,field_t,type_real,llm)70 CALL allocate_field(f_massflxc,field_u,type_real,llm)71 CALL allocate_field(f_massflx,field_u,type_real,llm)72 CALL allocate_field(f_uc,field_u,type_real,llm)73 first = .FALSE.74 END IF75 76 DO ind=1,ndomain77 CALL swap_dimensions(ind)78 CALL swap_geometry(ind)79 rhodz=f_rhodz(ind)80 massflx=f_massflx(ind)81 ps=f_ps(ind)82 u=f_u(ind)83 84 DO l = 1, llm85 DO j=jj_begin-1,jj_end+186 DO i=ii_begin-1,ii_end+187 ij=(j-1)*iim+i88 rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g89 ENDDO90 ENDDO91 ENDDO92 93 DO l = 1, llm94 DO j=jj_begin-1,jj_end+195 DO i=ii_begin-1,ii_end+196 ij=(j-1)*iim+i97 massflx(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)98 massflx(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)99 massflx(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)100 ENDDO101 ENDDO102 ENDDO103 ENDDO104 105 IF ( iadvtr == 0 ) THEN106 DO ind=1,ndomain107 CALL swap_dimensions(ind)108 CALL swap_geometry(ind)109 rhodz=f_rhodz(ind)110 rhodzm1 = f_rhodzm1(ind)111 massflxc = f_massflxc(ind) ! accumulated mass fluxes112 uc = f_uc(ind) ! time-integrated normal velocity to compute back-trajectory (Miura)113 rhodzm1 = rhodz114 massflxc = 0.0115 uc = 0.0116 END DO117 CALL transfert_request(f_rhodzm1,req_i1) !---->118 CALL transfert_request(f_massflxc,req_e1) !---->119 CALL transfert_request(f_massflxc,req_e1) !------>120 CALL transfert_request(f_uc,req_e1) !---->121 CALL transfert_request(f_uc,req_e1)122 END IF123 124 iadvtr = iadvtr + 1125 126 DO ind=1,ndomain127 CALL swap_dimensions(ind)128 CALL swap_geometry(ind)129 massflx = f_massflx(ind)130 massflxc = f_massflxc(ind)131 uc = f_uc(ind)132 u = f_u(ind)133 massflxc = massflxc + massflx134 uc = uc + u135 END DO136 137 IF ( iadvtr == iapp_tracvl ) THEN138 IF (is_mpi_root) PRINT *, 'Advection scheme'139 bigt = dt*iapp_tracvl140 DO ind=1,ndomain141 CALL swap_dimensions(ind)142 CALL swap_geometry(ind)143 uc = f_uc(ind)144 uc = uc/real(iapp_tracvl)145 massflxc = f_massflxc(ind)146 massflxc = massflxc*dt147 ! now massflx is time-integrated148 END DO149 150 CALL vlsplt(f_q,f_rhodzm1,f_massflxc,2.0,f_uc,bigt)151 iadvtr = 0152 END IF153 END SUBROUTINE advect_tracer154 155 SUBROUTINE vlsplt(f_q,f_rhodz,f_massflx,pente_max,f_u,bigt)156 58 USE field_mod 157 59 USE domain_mod … … 161 63 USE metric 162 64 USE advect_mod 65 USE advect_mod 66 USE disvert_mod 67 USE mpipara 163 68 IMPLICIT NONE 164 69 165 TYPE(t_field),POINTER :: f_q(:) 166 TYPE(t_field),POINTER :: f_u(:) 167 TYPE(t_field),POINTER :: f_rhodz(:) 168 TYPE(t_field),POINTER :: f_massflx(:) 169 170 TYPE(t_field),POINTER,SAVE :: f_wg(:) 171 TYPE(t_field),POINTER,SAVE :: f_zm(:) 172 TYPE(t_field),POINTER,SAVE :: f_zq(:) 173 174 REAL(rstd)::bigt 175 REAL(rstd),POINTER :: q(:,:,:) 176 REAL(rstd),POINTER :: u(:,:) 177 REAL(rstd),POINTER :: rhodz(:,:) 178 REAL(rstd),POINTER :: massflx(:,:) 179 REAL(rstd),POINTER :: wg(:,:) 180 REAL(rstd),POINTER :: zq(:,:,:) 181 REAL(rstd),POINTER :: zm(:,:) 182 REAL(rstd),POINTER :: tangent(:,:) 183 REAL(rstd),POINTER :: normal(:,:) 184 REAL(rstd),POINTER :: gradq3d(:,:,:) 185 186 REAL(rstd)::pente_max 187 LOGICAL,SAVE::first = .true. 188 INTEGER :: i,ij,l,j,ind,k 189 REAL(rstd) :: zzpbar, zzw 190 REAL::qvmax,qvmin 191 192 IF ( first ) THEN 193 CALL allocate_field(f_wg,field_t,type_real,llm) 194 CALL allocate_field(f_zm,field_t,type_real,llm) 195 CALL allocate_field(f_zq,field_t,type_real,llm,nqtot) 196 first = .FALSE. 197 END IF 70 TYPE(t_field),POINTER :: f_hfluxt(:) ! time-integrated horizontal mass flux 71 TYPE(t_field),POINTER :: f_wfluxt(:) ! time-integrated vertical mass flux 72 TYPE(t_field),POINTER :: f_u(:) ! velocity (for back-trajectories) 73 TYPE(t_field),POINTER :: f_q(:) ! tracer 74 TYPE(t_field),POINTER :: f_rhodz(:) ! mass field at beginning of macro time step 75 76 REAL(rstd),POINTER :: q(:,:,:), normal(:,:,:), tangent(:,:,:), gradq3d(:,:,:) 77 REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:) 78 REAL(rstd),POINTER :: rhodz(:,:), u(:,:) 79 INTEGER :: ind 80 81 CALL transfert_request(f_u,req_e1) 82 CALL transfert_request(f_u,req_e1) 83 CALL transfert_request(f_q,req_i1) 84 85 IF (is_mpi_root) PRINT *, 'Advection scheme' 198 86 199 87 DO ind=1,ndomain 200 88 CALL swap_dimensions(ind) 201 89 CALL swap_geometry(ind) 202 q=f_q(ind) 203 rhodz=f_rhodz(ind) 204 zq=f_zq(ind) 205 zm=f_zm(ind) 206 zm = rhodz ; zq = q 207 wg = f_wg(ind) 208 wg = 0.0 209 massflx=f_massflx(ind) 210 CALL advtrac(massflx,wg) ! compute vertical mass fluxes 211 END DO 212 213 DO ind=1,ndomain 214 CALL swap_dimensions(ind) 215 CALL swap_geometry(ind) 216 zq=f_zq(ind) 217 zm=f_zm(ind) 218 wg=f_wg(ind) 219 wg=wg*0.5 220 DO k = 1, nqtot 221 CALL vlz(zq(:,:,k),2.,zm,wg) 222 END DO 223 END DO 224 225 DO ind=1,ndomain 226 CALL swap_dimensions(ind) 227 CALL swap_geometry(ind) 228 zq=f_zq(ind) 229 zq = f_zq(ind) 230 zm = f_zm(ind) 231 massflx =f_massflx(ind) 232 u = f_u(ind) 90 q = f_q(ind) 91 rhodz = f_rhodz(ind) 92 hfluxt = f_hfluxt(ind) 93 wfluxt = f_wfluxt(ind) 94 normal = f_normal(ind) 233 95 tangent = f_tangent(ind) 234 normal = f_normal(ind)235 96 gradq3d = f_gradq3d(ind) 236 237 DO k = 1,nqtot 238 CALL compute_gradq3d(zq(:,:,k),gradq3d) 239 CALL compute_advect_horiz(tangent,normal,zq(:,:,k),gradq3d,zm,u,massflx,bigt) 240 END DO 241 END DO 242 243 DO ind=1,ndomain 244 CALL swap_dimensions(ind) 245 CALL swap_geometry(ind) 246 q = f_q(ind) 247 zq = f_zq(ind) 248 zm = f_zm(ind) 249 wg = f_wg(ind) 250 DO k = 1,nqtot 251 CALL vlz(zq(:,:,k),2.,zm,wg) 252 END DO 253 q = zq 254 END DO 255 256 END SUBROUTINE vlsplt 257 258 !============================================================================== 259 SUBROUTINE advtrac(massflx,wgg) 97 CALL compute_adv_tracer(tangent,normal,hfluxt,wfluxt,u, q,rhodz,gradq3d) 98 END DO 99 100 END SUBROUTINE advect_tracer 101 102 SUBROUTINE compute_adv_tracer(tangent,normal,hfluxt,wfluxt,u, q,rhodz,gradq3d) 103 USE icosa 104 USE field_mod 260 105 USE domain_mod 261 106 USE dimensions … … 263 108 USE geometry 264 109 USE metric 110 USE advect_mod 111 USE advect_mod 265 112 USE disvert_mod 113 USE mpipara 114 REAL(rstd), INTENT(IN) :: tangent(:,:,:), normal(:,:,:) 115 REAL(rstd), INTENT(IN) :: hfluxt(:,:), wfluxt(:,:), u(:,:) 116 REAL(rstd), INTENT(INOUT) :: rhodz(:,:) 117 REAL(rstd), INTENT(INOUT) :: q(:,:,:) 118 REAL(rstd), INTENT(OUT) :: gradq3d(:,:,:) 119 INTEGER :: k 120 ! for mass/tracer consistency each transport step also updates the mass field rhodz 121 ! mass is updated after all tracers have been updated, when k==nqtot 122 123 ! 1/2 vertical transport 124 DO k = 1, nqtot 125 CALL vlz(k==nqtot, 0.5,wfluxt,rhodz,q(:,:,k)) 126 END DO 127 128 ! horizontal transport 129 DO k = 1,nqtot 130 CALL compute_gradq3d(q(:,:,k),gradq3d) 131 CALL compute_advect_horiz(k==nqtot, tangent,normal,q(:,:,k),gradq3d,rhodz,u,hfluxt,dt*itau_adv) 132 END DO 133 134 ! 1/2 vertical transport 135 DO k = 1,nqtot 136 CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k)) 137 END DO 138 END SUBROUTINE compute_adv_tracer 139 140 SUBROUTINE vlz(update_mass, fac,wfluxt,mass, q) 141 ! 142 ! Auteurs: P.Le Van, F.Hourdin, F.Forget, T. Dubos 143 ! 144 ! ******************************************************************** 145 ! Update tracers using vertical mass flux only 146 ! Van Leer scheme with minmod limiter 147 ! wfluxt >0 for upward transport 148 ! ******************************************************************** 149 USE icosa 266 150 IMPLICIT NONE 267 REAL(rstd),INTENT(IN) :: massflx(iim*3*jjm,llm) 268 REAL(rstd),INTENT(OUT) :: wgg(iim*jjm,llm) 269 270 INTEGER :: i,j,ij,l 271 REAL(rstd) :: convm(iim*jjm,llm) 272 273 DO l = 1, llm 274 DO j=jj_begin,jj_end 275 DO i=ii_begin,ii_end 276 ij=(j-1)*iim+i 277 ! divergence of horizontal flux 278 convm(ij,l)= 1/(Ai(ij))*(ne(ij,right)*massflx(ij+u_right,l) + & 279 ne(ij,rup)*massflx(ij+u_rup,l) + & 280 ne(ij,lup)*massflx(ij+u_lup,l) + & 281 ne(ij,left)*massflx(ij+u_left,l) + & 282 ne(ij,ldown)*massflx(ij+u_ldown,l) + & 283 ne(ij,rdown)*massflx(ij+u_rdown,l)) 151 LOGICAL, INTENT(IN) :: update_mass 152 REAL(rstd), INTENT(IN) :: fac, wfluxt(iim*jjm,llm+1) ! vertical mass flux 153 REAL(rstd), INTENT(INOUT) :: mass(iim*jjm,llm) 154 REAL(rstd), INTENT(INOUT) :: q(iim*jjm,llm) 155 156 REAL(rstd) :: dq(iim*jjm,llm), & ! increase of q 157 dzqw(iim*jjm,llm), & ! vertical finite difference of q 158 adzqw(iim*jjm,llm), & ! abs(dzqw) 159 dzq(iim*jjm,llm), & ! limited slope of q 160 wq(iim*jjm,llm+1) ! time-integrated flux of q 161 REAL(rstd) :: dzqmax, newmass, sigw, qq, w 162 INTEGER :: i,ij,l,j 163 164 ! finite difference of q 165 DO l=2,llm 166 DO j=jj_begin,jj_end 167 DO i=ii_begin,ii_end 168 ij=(j-1)*iim+i 169 dzqw(ij,l)=q(ij,l)-q(ij,l-1) 170 adzqw(ij,l)=abs(dzqw(ij,l)) 284 171 ENDDO 285 172 ENDDO 286 173 ENDDO 287 174 288 ! accumulate divergence from top of model 289 DO l = llm-1, 1, -1 290 DO j=jj_begin,jj_end 291 DO i=ii_begin,ii_end 292 ij=(j-1)*iim+i 293 convm(ij,l) = convm(ij,l) + convm(ij,l+1) 175 ! minmod-limited slope of q 176 ! dzq = slope*dz, i.e. the reconstructed q varies by dzq inside level l 177 DO l=2,llm-1 178 DO j=jj_begin,jj_end 179 DO i=ii_begin,ii_end 180 ij=(j-1)*iim+i 181 IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN 182 dzq(ij,l) = 0.5*( dzqw(ij,l)+dzqw(ij,l+1) ) 183 dzqmax = pente_max * min( adzqw(ij,l),adzqw(ij,l+1) ) 184 dzq(ij,l) = sign( min(abs(dzq(ij,l)),dzqmax) , dzq(ij,l) ) ! NB : sign(a,b)=a*sign(b) 185 ELSE 186 dzq(ij,l)=0. 187 ENDIF 294 188 ENDDO 295 189 ENDDO 296 190 ENDDO 297 191 298 !!! Compute vertical velocity 299 DO l = 1,llm-1 300 DO j=jj_begin,jj_end 301 DO i=ii_begin,ii_end 302 ij=(j-1)*iim+i 303 wgg( ij, l+1 ) = (convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )) 304 ENDDO 305 ENDDO 306 ENDDO 307 192 ! 0 slope in top and bottom layers 308 193 DO j=jj_begin,jj_end 309 194 DO i=ii_begin,ii_end 310 195 ij=(j-1)*iim+i 311 wgg(ij,1) = 0. 196 dzq(ij,1)=0. 197 dzq(ij,llm)=0. 312 198 ENDDO 313 199 ENDDO 314 END SUBROUTINE advtrac 315 316 SUBROUTINE vlz(q,pente_max,masse,wgg) 317 !c 318 !c Auteurs: P.Le Van, F.Hourdin, F.Forget 319 !c 320 !c ******************************************************************** 321 !c Shema d'advection " pseudo amont " . 322 !c ******************************************************************** 323 USE icosa 324 IMPLICIT NONE 325 !c 326 !c Arguments: 327 !c ---------- 328 REAL masse(iim*jjm,llm),pente_max 329 REAL q(iim*jjm,llm) 330 REAL wgg(iim*jjm,llm),w(iim*jjm,llm+1) 331 REAL dq(iim*jjm,llm) 332 INTEGER i,ij,l,j,ii 333 !c 334 REAL wq(iim*jjm,llm+1),newmasse 335 336 REAL dzq(iim*jjm,llm),dzqw(iim*jjm,llm),adzqw(iim*jjm,llm),dzqmax 337 REAL sigw 338 339 REAL SSUM 340 341 342 w(:,1:llm) = -wgg(:,:) ! w>0 for downward transport 343 w(:,llm+1) = 0.0 344 345 !c On oriente tout dans le sens de la pression c'est a dire dans le 346 !c sens de W 347 348 DO l=2,llm 349 DO j=jj_begin,jj_end 350 DO i=ii_begin,ii_end 351 ij=(j-1)*iim+i 352 dzqw(ij,l)=q(ij,l-1)-q(ij,l) 353 adzqw(ij,l)=abs(dzqw(ij,l)) 354 ENDDO 355 ENDDO 356 ENDDO 357 358 DO l=2,llm-1 359 DO j=jj_begin,jj_end 360 DO i=ii_begin,ii_end 361 ij=(j-1)*iim+i 362 IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN 363 dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1)) 200 201 ! sigw = fraction of mass that leaves level l/l+1 202 ! then amount of q leaving level l/l+1 = wq = w * qq 203 DO l = 1,llm-1 204 DO j=jj_begin,jj_end 205 DO i=ii_begin,ii_end 206 ij=(j-1)*iim+i 207 IF(wfluxt(ij,l+1).gt.0.) THEN 208 ! downward transport, sigw<0 209 ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0 210 w = fac*wfluxt(ij,l+1) 211 sigw = w/mass(ij,l+1) 212 qq = q(ij,l+1) - 0.5*(1.+sigw)*dzq(ij,l+1) 213 wq(ij,l+1) = w*qq 364 214 ELSE 365 dzq(ij,l)=0. 215 ! upward transport, sigw>0 216 ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0 217 w = fac*wfluxt(ij,l) 218 sigw = w/mass(ij,l) 219 qq = q(ij,l)+0.5*(1.-sigw)*dzq(ij,l) 220 wq(ij,l+1) = w*qq 366 221 ENDIF 367 dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1)) 368 dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l)) 369 ENDDO 370 ENDDO 371 ENDDO 372 373 DO l=2,llm-1 374 DO j=jj_begin,jj_end 375 DO i=ii_begin,ii_end 376 ij=(j-1)*iim+i 377 dzq(ij,1)=0. 378 dzq(ij,llm)=0. 379 ENDDO 380 ENDDO 381 ENDDO 382 !c --------------------------------------------------------------- 383 !c .... calcul des termes d'advection verticale ....... 384 !c --------------------------------------------------------------- 385 386 !c calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq 387 388 DO l = 1,llm-1 389 DO j=jj_begin,jj_end 390 DO i=ii_begin,ii_end 391 ij=(j-1)*iim+i 392 IF(w(ij,l+1).gt.0.) THEN 393 ! upwind only if downward transport 394 sigw=w(ij,l+1)/masse(ij,l+1) 395 wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1)) 396 ELSE 397 ! upwind only if upward transport 398 sigw=w(ij,l+1)/masse(ij,l) 399 wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l)) 400 ENDIF 401 ENDDO 402 ENDDO 403 END DO 404 222 ENDDO 223 ENDDO 224 END DO 225 ! wq = 0 at top and bottom 405 226 DO j=jj_begin,jj_end 406 227 DO i=ii_begin,ii_end … … 411 232 END DO 412 233 234 ! update q, mass is updated only after all q's have been updated 413 235 DO l=1,llm 414 236 DO j=jj_begin,jj_end 415 237 DO i=ii_begin,ii_end 416 238 ij=(j-1)*iim+i 417 ! masse -= dw/dz but w>0 <=> downward 418 newmasse=masse(ij,l)+(w(ij,l+1)-w(ij,l)) 419 ! dq(ij,l) = (wq(ij,l+1)-wq(ij,l)) !================>>>> 420 q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))/newmasse 421 ! masse(ij,l)=newmasse 422 ENDDO 423 ENDDO 424 END DO 425 RETURN 239 newmass = mass(ij,l) + fac*(wfluxt(ij,l)-wfluxt(ij,l+1)) 240 q(ij,l) = ( q(ij,l)*mass(ij,l) + wq(ij,l)-wq(ij,l+1) ) / newmass 241 IF(update_mass) mass(ij,l)=newmass 242 ENDDO 243 ENDDO 244 END DO 245 426 246 END SUBROUTINE vlz 427 247
Note: See TracChangeset
for help on using the changeset viewer.