Changeset 181
- Timestamp:
- 12/15/13 01:54:31 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/caldyn_gcm.f90
r179 r181 426 426 427 427 IF(caldyn_eta==eta_mass) THEN 428 428 429 429 !!! Compute exner function and geopotential 430 430 DO l = 1,llm 431 431 !$OMP DO SCHEDULE(STATIC) 432 !DIR$ SIMD433 DO ij=ij_begin_ext,ij_end_ext434 p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later435 ! p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j))436 exner_ik = cpp * (p_ik/preff) ** kappa437 pk(ij,l) = exner_ik438 ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz439 geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik432 !DIR$ SIMD 433 DO ij=ij_begin_ext,ij_end_ext 434 p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later 435 ! p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j)) 436 exner_ik = cpp * (p_ik/preff) ** kappa 437 pk(ij,l) = exner_ik 438 ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 439 geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 440 440 ENDDO 441 441 ENDDO … … 447 447 ! Notice that the computation below should work also when caldyn_eta=eta_mass 448 448 449 ! uppermost layer 450 !DIR$ SIMD 451 DO ij=ij_begin_ext,ij_end_ext 452 pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm) 453 END DO 454 ! other layers 455 DO l = llm-1, 1, -1 456 !$OMP DO SCHEDULE(STATIC) 457 !DIR$ SIMD 449 IF(boussinesq) THEN 450 ! use theta*rhodz in hydrostatic balance 451 ! uppermost layer 452 !DIR$ SIMD 458 453 DO ij=ij_begin_ext,ij_end_ext 459 pk(ij,l ) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1))454 pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm) 460 455 END DO 461 END DO 462 ! surface pressure (for diagnostics) 463 DO ij=ij_begin_ext,ij_end_ext 464 ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1) 465 END DO 466 467 IF(boussinesq) THEN ! specific volume 1 = dphi/g/rhodz 456 ! other layers 457 DO l = llm-1, 1, -1 458 !$OMP DO SCHEDULE(STATIC) 459 !DIR$ SIMD 460 DO ij=ij_begin_ext,ij_end_ext 461 pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1)) 462 END DO 463 END DO 464 ! surface pressure (for diagnostics) 465 DO ij=ij_begin_ext,ij_end_ext 466 ps(ij) = pk(ij,1) + (.5*g)*theta(ij,1)*rhodz(ij,1) 467 END DO 468 ! now pk contains the Lagrange multiplier (pressure) 469 470 ! specific volume 1 = dphi/g/rhodz 468 471 DO l = 1,llm 469 472 !$OMP DO SCHEDULE(STATIC) 470 !DIR$ SIMD473 !DIR$ SIMD 471 474 DO ij=ij_begin_ext,ij_end_ext 472 geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)475 geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l) 473 476 ENDDO 474 477 ENDDO 475 ELSE ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 478 ELSE ! non-Boussinesq 479 ! uppermost layer 480 !DIR$ SIMD 481 DO ij=ij_begin_ext,ij_end_ext 482 pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm) 483 END DO 484 ! other layers 485 DO l = llm-1, 1, -1 486 !$OMP DO SCHEDULE(STATIC) 487 !DIR$ SIMD 488 DO ij=ij_begin_ext,ij_end_ext 489 pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1)) 490 END DO 491 END DO 492 ! surface pressure (for diagnostics) 493 DO ij=ij_begin_ext,ij_end_ext 494 ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1) 495 END DO 496 497 498 ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 476 499 DO l = 1,llm 477 500 !$OMP DO SCHEDULE(STATIC) 478 !DIR$ SIMD479 501 !DIR$ SIMD 502 DO ij=ij_begin_ext,ij_end_ext 480 503 p_ik = pk(ij,l) 481 504 exner_ik = cpp * (p_ik/preff) ** kappa … … 666 689 !DIR$ SIMD 667 690 DO ij=ij_begin,ij_end 668 691 ! until this point pk contains the Lagrange multiplier (pressure) 669 692 berni(ij,l) = pk(ij,l) & 670 693 + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 + & … … 674 697 le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 + & 675 698 le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 699 ! from now on pk contains the vertically-averaged geopotential 676 700 pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) 677 701 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.