Changeset 183 for codes/icosagcm/trunk
- Timestamp:
- 12/16/13 00:48:24 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/caldyn_gcm.f90
r182 r183 221 221 pk = f_pk(ind) 222 222 geopot = f_geopot(ind) 223 CALL compute_geopot(ps,mass,theta, u,pk,geopot)223 CALL compute_geopot(ps,mass,theta, pk,geopot) 224 224 hflux=f_hflux(ind) 225 225 convm = f_dmass(ind) … … 249 249 pk = f_pk(ind) 250 250 geopot = f_geopot(ind) 251 CALL compute_geopot(ps,mass,theta, u,pk,geopot)251 CALL compute_geopot(ps,mass,theta, pk,geopot) 252 252 hflux=f_hflux(ind) 253 253 convm = f_dmass(ind) … … 405 405 END SUBROUTINE compute_pvort 406 406 407 SUBROUTINE compute_geopot(ps,rhodz,theta, u,pk,geopot)407 SUBROUTINE compute_geopot(ps,rhodz,theta, pk,geopot) 408 408 USE icosa 409 409 USE disvert_mod … … 415 415 REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 416 416 REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm) ! potential temperature 417 REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm) ! prognostic velocity418 417 REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function 419 418 REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential 420 419 421 REAL(rstd) :: urel(3*iim*jjm,llm) ! relative velocity422 REAL(rstd) :: cor_NT(iim*jjm,llm) ! NT coriolis force u.(du/dPhi)423 420 INTEGER :: i,j,ij,l 424 421 REAL(rstd) :: p_ik, exner_ik … … 448 445 ! Notice that the computation below should work also when caldyn_eta=eta_mass 449 446 450 IF(boussinesq) THEN 451 ! use theta*rhodz in hydrostatic balance 452 ! uppermost layer 453 !DIR$ SIMD 454 DO ij=ij_begin_ext,ij_end_ext 455 pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm) 456 END DO 457 ! other layers 458 DO l = llm-1, 1, -1 459 !$OMP DO SCHEDULE(STATIC) 460 !DIR$ SIMD 461 DO ij=ij_begin_ext,ij_end_ext 462 pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1)) 463 END DO 464 END DO 465 ! surface pressure (for diagnostics) 466 DO ij=ij_begin_ext,ij_end_ext 467 ps(ij) = pk(ij,1) + (.5*g)*theta(ij,1)*rhodz(ij,1) 468 END DO 469 ! now pk contains the Lagrange multiplier (pressure) 470 447 IF(boussinesq) THEN ! compute only geopotential : pressure pk will be computed in compute_caldyn_horiz 471 448 ! specific volume 1 = dphi/g/rhodz 472 449 DO l = 1,llm … … 477 454 ENDDO 478 455 ENDDO 479 ELSE ! non-Boussinesq 456 ELSE ! non-Boussinesq, compute geopotential and Exner pressure 480 457 ! uppermost layer 481 458 !DIR$ SIMD … … 496 473 END DO 497 474 498 499 475 ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 500 476 DO l = 1,llm … … 523 499 USE omp_para 524 500 IMPLICIT NONE 525 REAL(rstd),INTENT(IN) :: u(iim*3*jjm,llm) 501 REAL(rstd),INTENT(IN) :: u(iim*3*jjm,llm) ! prognostic "velocity" 526 502 REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 527 503 REAL(rstd),INTENT(IN) :: qu(iim*3*jjm,llm) … … 535 511 REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm) 536 512 513 REAL(rstd) :: cor_NT(iim*jjm,llm) ! NT coriolis force u.(du/dPhi) 537 514 REAL(rstd) :: urel(3*iim*jjm,llm) ! relative velocity 538 515 REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux … … 686 663 !!! Compute bernouilli term = Kinetic Energy + geopotential 687 664 IF(boussinesq) THEN 665 ! first use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure) 666 ! uppermost layer 667 !DIR$ SIMD 668 DO ij=ij_begin_ext,ij_end_ext 669 pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm) 670 END DO 671 ! other layers 672 DO l = llm-1, 1, -1 673 !$OMP DO SCHEDULE(STATIC) 674 !DIR$ SIMD 675 DO ij=ij_begin_ext,ij_end_ext 676 pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1)) 677 END DO 678 END DO 679 ! surface pressure (for diagnostics) FIXME 680 ! DO ij=ij_begin_ext,ij_end_ext 681 ! ps(ij) = pk(ij,1) + (.5*g)*theta(ij,1)*rhodz(ij,1) 682 ! END DO 683 ! now pk contains the Lagrange multiplier (pressure) 688 684 689 685 DO l=ll_begin,ll_end 690 686 !DIR$ SIMD 691 687 DO ij=ij_begin,ij_end 692 ! until this point pk contains the Lagrange multiplier (pressure)693 688 berni(ij,l) = pk(ij,l) & 694 689 + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 + & … … 703 698 ENDDO 704 699 705 ELSE 700 ELSE ! compressible 706 701 707 702 DO l=ll_begin,ll_end
Note: See TracChangeset
for help on using the changeset viewer.