Changeset 134 for codes/icosagcm/trunk
- Timestamp:
- 02/14/13 18:26:54 (11 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/caldyn.f90
r129 r134 20 20 CASE('gcm') 21 21 CALL init_caldyn_gcm 22 CASE('adv')23 CALL init_caldyn_adv22 ! CASE('adv') 23 ! CALL init_caldyn_adv 24 24 CASE DEFAULT 25 25 PRINT*,'Bad selector for variable caldyn : <', TRIM(caldyn_type),'> options are <gcm>, <adv>' … … 30 30 END SUBROUTINE init_caldyn 31 31 32 SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 32 SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, & 33 f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) 33 34 USE icosa 34 35 USE caldyn_gcm_mod, ONLY : caldyn_gcm=>caldyn … … 41 42 TYPE(t_field),POINTER :: f_u(:) 42 43 TYPE(t_field),POINTER :: f_q(:) 44 TYPE(t_field),POINTER :: f_hflux(:) 45 TYPE(t_field),POINTER :: f_wflux(:) 43 46 TYPE(t_field),POINTER :: f_dps(:) 44 47 TYPE(t_field),POINTER :: f_dtheta_rhodz(:) … … 47 50 SELECT CASE (TRIM(caldyn_type)) 48 51 CASE('gcm') 49 CALL caldyn_gcm(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 50 CASE('adv') 51 CALL caldyn_adv(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 52 CALL caldyn_gcm(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, & 53 f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) 54 ! CASE('adv') 55 ! CALL caldyn_adv(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 52 56 END SELECT 53 57 -
codes/icosagcm/trunk/src/caldyn_gcm.f90
r133 r134 35 35 STOP 36 36 END SELECT 37 37 IF (is_mpi_root) PRINT *, 'caldyn_conserv=',def 38 38 39 39 def='direct' … … 84 84 END SUBROUTINE allocate_caldyn 85 85 86 SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du) 86 SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, & 87 f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) 87 88 USE icosa 88 89 USE vorticity_mod … … 97 98 TYPE(t_field),POINTER :: f_u(:) 98 99 TYPE(t_field),POINTER :: f_q(:) 100 TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:) 99 101 TYPE(t_field),POINTER :: f_dps(:) 100 102 TYPE(t_field),POINTER :: f_dtheta_rhodz(:) … … 103 105 REAL(rstd),POINTER :: phis(:), ps(:), dps(:) 104 106 REAL(rstd),POINTER :: theta_rhodz(:,:), dtheta_rhodz(:,:) 105 REAL(rstd),POINTER :: u(:,:), du(:,:) 107 REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:) 106 108 REAL(rstd),POINTER :: p(:,:), rhodz(:,:), qu(:,:) 107 109 INTEGER :: ind,ij … … 133 135 CALL swap_geometry(ind) 134 136 phis=f_phis(ind) 137 hflux=f_hflux(ind) 138 wflux=f_wflux(ind) 135 139 ps=f_ps(ind) 136 140 dps=f_dps(ind) … … 144 148 out_u=f_out_u(ind) 145 149 !$OMP PARALLEL DEFAULT(SHARED) 146 CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du) 150 CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, & 151 hflux, wflux, dps, dtheta_rhodz, du) 147 152 !$OMP END PARALLEL 148 153 ENDDO … … 155 160 ps=f_ps(ind) 156 161 dps=f_dps(ind) 162 hflux=f_hflux(ind) 163 wflux=f_wflux(ind) 157 164 theta_rhodz=f_theta_rhodz(ind) 158 165 dtheta_rhodz=f_dtheta_rhodz(ind) … … 165 172 !$OMP PARALLEL DEFAULT(SHARED) 166 173 CALL compute_pvort(ps, u, p,rhodz,qu) 167 CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du) 174 CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, & 175 hflux, wflux, dps, dtheta_rhodz, du) 168 176 !$OMP END PARALLEL 169 177 ENDDO … … 276 284 END SUBROUTINE compute_pvort 277 285 278 SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du)286 SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du) 279 287 USE icosa 280 288 USE disvert_mod … … 289 297 REAL(rstd),INTENT(IN) :: qu(iim*3*jjm,llm) 290 298 291 REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm) 299 REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm), hflux(iim*3*jjm,llm) 292 300 REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm) 293 301 REAL(rstd),INTENT(OUT) :: dps(iim*jjm) 302 REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux 294 303 295 304 INTEGER :: i,j,ij,l … … 301 310 REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:), pks(:) ! Exner function 302 311 REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:), beta(:,:) 303 REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:) 304 REAL(rstd),ALLOCATABLE,SAVE :: F e(:,:), Ftheta(:,:) ! mass flux,theta flux312 REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:) ! geopotential 313 REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux 305 314 REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:) ! mass flux convergence 306 REAL(rstd),ALLOCATABLE,SAVE :: w(:,:) ! vertical mass flux307 315 REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:) ! Bernouilli function 308 316 … … 319 327 ALLOCATE(beta(iim*jjm,llm)) 320 328 ALLOCATE(phi(iim*jjm,llm)) ! geopotential 321 ALLOCATE(Fe(3*iim*jjm,llm)) ! mass flux322 329 ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux 323 330 ALLOCATE(convm(iim*jjm,llm)) ! mass flux convergence 324 ALLOCATE(w(iim*jjm,llm)) ! vertical velocity325 331 ALLOCATE(berni(iim*jjm,llm)) ! bernouilli term 326 332 … … 342 348 DO i=ii_begin-1,ii_end+1 343 349 ij=(j-1)*iim+i 344 Fe(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)345 Fe(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)346 Fe(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)347 348 Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))* Fe(ij+u_right,l)349 Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))* Fe(ij+u_lup,l)350 Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))* Fe(ij+u_ldown,l)350 hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) 351 hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) 352 hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown) 353 354 Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*hflux(ij+u_right,l) 355 Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l) 356 Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*hflux(ij+u_ldown,l) 351 357 ENDDO 352 358 ENDDO … … 358 364 ij=(j-1)*iim+i 359 365 ! convm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 360 convm(ij,l)= 1./Ai(ij)*(ne(ij,right)* Fe(ij+u_right,l) + &361 ne(ij,rup)* Fe(ij+u_rup,l) + &362 ne(ij,lup)* Fe(ij+u_lup,l) + &363 ne(ij,left)* Fe(ij+u_left,l) + &364 ne(ij,ldown)* Fe(ij+u_ldown,l) + &365 ne(ij,rdown)* Fe(ij+u_rdown,l))366 convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*hflux(ij+u_right,l) + & 367 ne(ij,rup)*hflux(ij+u_rup,l) + & 368 ne(ij,lup)*hflux(ij+u_lup,l) + & 369 ne(ij,left)*hflux(ij+u_left,l) + & 370 ne(ij,ldown)*hflux(ij+u_ldown,l) + & 371 ne(ij,rdown)*hflux(ij+u_rdown,l)) 366 372 367 373 ! signe ? attention d (rho theta dz) … … 396 402 ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt 397 403 ! => w>0 for upward transport 398 w ( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )404 wflux( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 ) 399 405 ENDDO 400 406 ENDDO … … 406 412 DO i=ii_begin,ii_end 407 413 ij=(j-1)*iim+i 408 w (ij,1) = 0.414 wflux(ij,1) = 0. 409 415 ! dps/dt = -int(div flux)dz 410 416 dps(ij)=-convm(ij,1) * g … … 423 429 ij=(j-1)*iim+i 424 430 425 uu = wee(ij+u_right,1,1)* Fe(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+ &426 wee(ij+u_right,2,1)* Fe(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+ &427 wee(ij+u_right,3,1)* Fe(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+ &428 wee(ij+u_right,4,1)* Fe(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+ &429 wee(ij+u_right,5,1)* Fe(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+ &430 wee(ij+u_right,1,2)* Fe(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+ &431 wee(ij+u_right,2,2)* Fe(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+ &432 wee(ij+u_right,3,2)* Fe(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+ &433 wee(ij+u_right,4,2)* Fe(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+ &434 wee(ij+u_right,5,2)* Fe(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))431 uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+ & 432 wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+ & 433 wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+ & 434 wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+ & 435 wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+ & 436 wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+ & 437 wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+ & 438 wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+ & 439 wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+ & 440 wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l)) 435 441 du(ij+u_right,l) = .5*uu/de(ij+u_right) 436 442 437 uu = wee(ij+u_lup,1,1)* Fe(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) + &438 wee(ij+u_lup,2,1)* Fe(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) + &439 wee(ij+u_lup,3,1)* Fe(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) + &440 wee(ij+u_lup,4,1)* Fe(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) + &441 wee(ij+u_lup,5,1)* Fe(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) + &442 wee(ij+u_lup,1,2)* Fe(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) + &443 wee(ij+u_lup,2,2)* Fe(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) + &444 wee(ij+u_lup,3,2)* Fe(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) + &445 wee(ij+u_lup,4,2)* Fe(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) + &446 wee(ij+u_lup,5,2)* Fe(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))443 uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) + & 444 wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) + & 445 wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) + & 446 wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) + & 447 wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) + & 448 wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) + & 449 wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) + & 450 wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) + & 451 wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) + & 452 wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l)) 447 453 du(ij+u_lup,l) = .5*uu/de(ij+u_lup) 448 454 449 455 450 uu = wee(ij+u_ldown,1,1)* Fe(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) + &451 wee(ij+u_ldown,2,1)* Fe(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) + &452 wee(ij+u_ldown,3,1)* Fe(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) + &453 wee(ij+u_ldown,4,1)* Fe(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) + &454 wee(ij+u_ldown,5,1)* Fe(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) + &455 wee(ij+u_ldown,1,2)* Fe(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) + &456 wee(ij+u_ldown,2,2)* Fe(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) + &457 wee(ij+u_ldown,3,2)* Fe(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) + &458 wee(ij+u_ldown,4,2)* Fe(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) + &459 wee(ij+u_ldown,5,2)* Fe(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))456 uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) + & 457 wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) + & 458 wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) + & 459 wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) + & 460 wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) + & 461 wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) + & 462 wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) + & 463 wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) + & 464 wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) + & 465 wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l)) 460 466 du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown) 461 467 … … 472 478 ij=(j-1)*iim+i 473 479 474 uu = wee(ij+u_right,1,1)* Fe(ij+u_rup,l)+ &475 wee(ij+u_right,2,1)* Fe(ij+u_lup,l)+ &476 wee(ij+u_right,3,1)* Fe(ij+u_left,l)+ &477 wee(ij+u_right,4,1)* Fe(ij+u_ldown,l)+ &478 wee(ij+u_right,5,1)* Fe(ij+u_rdown,l)+ &479 wee(ij+u_right,1,2)* Fe(ij+t_right+u_ldown,l)+ &480 wee(ij+u_right,2,2)* Fe(ij+t_right+u_rdown,l)+ &481 wee(ij+u_right,3,2)* Fe(ij+t_right+u_right,l)+ &482 wee(ij+u_right,4,2)* Fe(ij+t_right+u_rup,l)+ &483 wee(ij+u_right,5,2)* Fe(ij+t_right+u_lup,l)480 uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+ & 481 wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+ & 482 wee(ij+u_right,3,1)*hflux(ij+u_left,l)+ & 483 wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+ & 484 wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+ & 485 wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+ & 486 wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+ & 487 wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+ & 488 wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+ & 489 wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l) 484 490 du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right) 485 491 486 492 487 uu = wee(ij+u_lup,1,1)* Fe(ij+u_left,l)+ &488 wee(ij+u_lup,2,1)* Fe(ij+u_ldown,l)+ &489 wee(ij+u_lup,3,1)* Fe(ij+u_rdown,l)+ &490 wee(ij+u_lup,4,1)* Fe(ij+u_right,l)+ &491 wee(ij+u_lup,5,1)* Fe(ij+u_rup,l)+ &492 wee(ij+u_lup,1,2)* Fe(ij+t_lup+u_right,l)+ &493 wee(ij+u_lup,2,2)* Fe(ij+t_lup+u_rup,l)+ &494 wee(ij+u_lup,3,2)* Fe(ij+t_lup+u_lup,l)+ &495 wee(ij+u_lup,4,2)* Fe(ij+t_lup+u_left,l)+ &496 wee(ij+u_lup,5,2)* Fe(ij+t_lup+u_ldown,l)493 uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+ & 494 wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+ & 495 wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+ & 496 wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+ & 497 wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+ & 498 wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+ & 499 wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+ & 500 wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+ & 501 wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+ & 502 wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l) 497 503 du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup) 498 504 499 uu = wee(ij+u_ldown,1,1)* Fe(ij+u_rdown,l)+ &500 wee(ij+u_ldown,2,1)* Fe(ij+u_right,l)+ &501 wee(ij+u_ldown,3,1)* Fe(ij+u_rup,l)+ &502 wee(ij+u_ldown,4,1)* Fe(ij+u_lup,l)+ &503 wee(ij+u_ldown,5,1)* Fe(ij+u_left,l)+ &504 wee(ij+u_ldown,1,2)* Fe(ij+t_ldown+u_lup,l)+ &505 wee(ij+u_ldown,2,2)* Fe(ij+t_ldown+u_left,l)+ &506 wee(ij+u_ldown,3,2)* Fe(ij+t_ldown+u_ldown,l)+ &507 wee(ij+u_ldown,4,2)* Fe(ij+t_ldown+u_rdown,l)+ &508 wee(ij+u_ldown,5,2)* Fe(ij+t_ldown+u_right,l)505 uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+ & 506 wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+ & 507 wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+ & 508 wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+ & 509 wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+ & 510 wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+ & 511 wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+ & 512 wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+ & 513 wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+ & 514 wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l) 509 515 du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown) 510 516 … … 606 612 ! ww>0 <=> upward transport 607 613 608 ww = 0.5 * w (ij,l+1) * (theta(ij,l) + theta(ij,l+1) )614 ww = 0.5 * wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1) ) 609 615 dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l ) - ww 610 616 dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1) + ww 611 617 612 ww = 0.5 * ( w (ij,l+1) + w(ij+t_right,l+1))618 ww = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_right,l+1)) 613 619 uu = u(ij+u_right,l+1) - u(ij+u_right,l) 614 620 du(ij+u_right, l ) = du(ij+u_right,l) - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))) 615 621 du(ij+u_right, l+1 ) = du(ij+u_right,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_right,l+1))) 616 622 617 ww = 0.5 * ( w (ij,l+1) + w(ij+t_lup,l+1))623 ww = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_lup,l+1)) 618 624 uu = u(ij+u_lup,l+1) - u(ij+u_lup,l) 619 625 du(ij+u_lup, l ) = du(ij+u_lup,l) - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))) 620 626 du(ij+u_lup, l+1 ) = du(ij+u_lup,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_lup,l+1))) 621 627 622 ww = 0.5 * ( w (ij,l+1) + w(ij+t_ldown,l+1))628 ww = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_ldown,l+1)) 623 629 uu = u(ij+u_ldown,l+1) - u(ij+u_ldown,l) 624 630 du(ij+u_ldown, l ) = du(ij+u_ldown,l) - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))) … … 631 637 !!$OMP BARRIER 632 638 !!$OMP MASTER 633 DEALLOCATE(theta) ! potential temperature634 DEALLOCATE(pk) ! Exner function639 DEALLOCATE(theta) ! potential temperature 640 DEALLOCATE(pk) ! Exner function 635 641 DEALLOCATE(pks) 636 642 DEALLOCATE(alpha) 637 643 DEALLOCATE(beta) 638 DEALLOCATE(phi) ! geopotential 639 DEALLOCATE(Fe) ! mass flux 640 DEALLOCATE(Ftheta) ! theta flux 641 DEALLOCATE(convm) ! mass flux convergence 642 DEALLOCATE(w) ! vertical velocity 644 DEALLOCATE(phi) ! geopotential 645 DEALLOCATE(Ftheta) ! theta flux 646 DEALLOCATE(convm) ! mass flux convergence 643 647 DEALLOCATE(berni) ! bernouilli term 644 648 !!$OMP END MASTER -
codes/icosagcm/trunk/src/timeloop_gcm.f90
r133 r134 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 time 48 49 ! INTEGER :: itaumax 49 50 ! REAL(rstd) ::write_period … … 78 79 CALL allocate_field(f_q,field_t,type_real,llm,nqtot) 79 80 CALL allocate_field(f_rhodz,field_t,type_real,llm) 81 ! Mass fluxes 80 82 CALL allocate_field(f_hflux,field_u,type_real,llm) ! instantaneous mass fluxes 81 83 CALL allocate_field(f_wflux,field_t,type_real,llm+1) ! computed by caldyn … … 138 140 CALL swap_dimensions(ind) 139 141 CALL swap_geometry(ind) 140 hfluxt=f_hfluxt(ind); hfluxt=0. ! clear accumulated fluxes141 wfluxt=f_hfluxt(ind); wfluxt=0.142 142 rhodz=f_rhodz(ind); ps=f_ps(ind) 143 143 CALL compute_rhodz(ps,rhodz) ! save rhodz for transport scheme before dynamics update ps 144 144 END DO 145 fluxt_zero=.FALSE. 145 146 146 147 DO it=0,itaumax … … 155 156 156 157 DO stage=1,nb_stage 157 ! CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &158 ! f_phis,f_ps,f_theta_rhodz,f_u, f_q, &159 ! f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)160 158 CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), & 161 159 f_phis,f_ps,f_theta_rhodz,f_u, f_q, & 162 f_ dps, f_dtheta_rhodz, f_du)160 f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) 163 161 SELECT CASE (scheme) 164 162 CASE(euler) … … 185 183 IF(MOD(it+1,itau_adv)==0) THEN 186 184 ! CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz) ! update q and rhodz after RK step 187 DO ind=1,ndomain 188 hfluxt=f_hfluxt(ind) 189 wfluxt=f_hfluxt(ind) 190 hfluxt=0. ! clear accumulated fluxes 191 wfluxt=0. 192 END DO 185 fluxt_zero=.TRUE. 193 186 END IF 194 187 … … 209 202 hflux=f_hflux(ind); hfluxt=f_hfluxt(ind) 210 203 wflux=f_hflux(ind); wfluxt=f_wfluxt(ind) 211 CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt )204 CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero) 212 205 END IF 213 206 u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) … … 242 235 hflux=f_hflux(ind); hfluxt=f_hfluxt(ind) 243 236 wflux=f_hflux(ind); wfluxt=f_wfluxt(ind) 244 CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt )237 CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero) 245 238 END IF 246 239 END DO … … 366 359 END SUBROUTINE compute_rhodz 367 360 368 SUBROUTINE accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,tau )361 SUBROUTINE accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,tau,fluxt_zero) 369 362 USE icosa 370 363 REAL(rstd), INTENT(IN) :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1) 371 364 REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1) 372 365 REAL(rstd), INTENT(IN) :: tau 373 hfluxt = hfluxt + tau*hflux 374 wfluxt = wfluxt + tau*wflux 366 LOGICAL, INTENT(INOUT) :: fluxt_zero 367 IF(fluxt_zero) THEN 368 fluxt_zero=.FALSE. 369 hfluxt = tau*hflux 370 wfluxt = tau*wflux 371 ELSE 372 hfluxt = hfluxt + tau*hflux 373 wfluxt = wfluxt + tau*wflux 374 END IF 375 375 END SUBROUTINE accumulate_fluxes 376 376
Note: See TracChangeset
for help on using the changeset viewer.