Changeset 580
- Timestamp:
- 10/13/17 15:55:07 (7 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 12 added
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/base/earth_const.f90
r553 r580 28 28 !$OMP THREADPRIVATE(boussinesq) 29 29 !$OMP THREADPRIVATE(hydrostatic) 30 LOGICAL :: dysl, dysl_geopot, dysl_pvort_only, dysl_caldyn_fast, dysl_caldyn_coriolis, dysl_slow_hydro 31 !$OMP THREADPRIVATE(dysl, dysl_geopot, dysl_pvort_only, dysl_caldyn_fast, dysl_caldyn_coriolis, dysl_slow_hydro) 30 32 31 33 CONTAINS -
codes/icosagcm/trunk/src/dynamics/caldyn.f90
r548 r580 2 2 USE icosa 3 3 PRIVATE 4 SAVE 4 5 CHARACTER(LEN=255),SAVE :: caldyn_type 5 6 !$OMP THREADPRIVATE(caldyn_type) 6 7 PUBLIC init_caldyn, caldyn, caldyn_BC 7 8 PUBLIC init_caldyn, caldyn, caldyn_BC, dysl 8 9 9 10 CONTAINS … … 28 29 END SELECT 29 30 30 31 31 END SUBROUTINE init_caldyn 32 32 -
codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90
r561 r580 22 22 REAL(rstd),POINTER :: planetvel(:) 23 23 24 #ifdef CPP_DYSL25 IF(is_master) PRINT *,'CPP_DYSL : Using macro-generated compute kernels'26 #endif27 28 24 hydrostatic=.TRUE. 29 25 CALL getin("hydrostatic",hydrostatic) 30 26 27 dysl=.NOT.hydrostatic ! dysl defaults to .FALSE. if hydrostatic, .TRUE. if NH 28 CALL getin("dysl",dysl) 29 dysl_geopot=dysl 30 CALL getin("dysl_geopot",dysl_geopot) 31 dysl_caldyn_fast=dysl 32 CALL getin("dysl_caldyn_fast",dysl_caldyn_fast) 33 dysl_slow_hydro=dysl 34 CALL getin("dysl_slow_hydro",dysl_slow_hydro) 35 dysl_pvort_only=dysl 36 CALL getin("dysl_pvort_only",dysl_pvort_only) 37 dysl_caldyn_coriolis=dysl 38 CALL getin("dysl_caldyn_coriolis",dysl_caldyn_coriolis) 39 31 40 def='energy' 32 41 CALL getin('caldyn_conserv',def) -
codes/icosagcm/trunk/src/dynamics/caldyn_hevi.f90
r561 r580 48 48 TYPE(t_field) :: f_dPhi_fast(:) 49 49 50 REAL(rstd),POINTER :: ps(:), dps(:) 50 REAL(rstd),POINTER :: ps(:), dps(:), phis(:) 51 51 REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:) 52 52 REAL(rstd),POINTER :: du(:,:), dW(:,:), dPhi(:,:), hflux(:,:), wflux(:,:) … … 115 115 CALL compute_geopot(mass,theta, ps,pk,geopot) 116 116 ELSE 117 phis = f_phis(ind) 117 118 W = f_W(ind) 118 119 dW = f_dW_fast(ind) … … 121 122 m_il = f_wil(ind) 122 123 pres = f_gradPhi2(ind) 123 CALL compute_caldyn_solver(tau, mass,theta,pk,geopot,W, m_il,pres, dPhi,dW,du) ! computes d(Phi,W,du)_fast and updates Phi,W124 CALL compute_caldyn_solver(tau,phis, mass,theta,pk,geopot,W, m_il,pres, dPhi,dW,du) ! computes d(Phi,W,du)_fast and updates Phi,W 124 125 END IF 125 126 u=f_u(ind) -
codes/icosagcm/trunk/src/dynamics/caldyn_kernels_base.F90
r561 r580 49 49 Rd = kappa*cpp 50 50 51 #ifdef CPP_DYSL 52 !#if 0 51 IF(dysl_geopot) THEN 53 52 #include "../kernels/compute_geopot.k90" 54 #else 55 53 ELSE 56 54 ! Pressure is computed first top-down (temporarily stored in pk) 57 55 ! Then Exner pressure and geopotential are computed bottom-up … … 164 162 END IF 165 163 166 #endif 164 END IF ! dysl 167 165 168 166 !ym flush geopot … … 314 312 CALL trace_start("compute_caldyn_vert_nh") 315 313 316 #ifdef CPP_DYSL 314 IF(dysl) THEN 317 315 !$OMP BARRIER 318 316 #include "../kernels/caldyn_vert_NH.k90" 319 317 !$OMP BARRIER 320 #else 318 ELSE 321 319 #define ETA_DOT(ij) eta_dot(ij,1) 322 320 #define WCOV(ij) wcov(ij,1) … … 365 363 #undef ETA_DOT 366 364 #undef WCOV 367 #endif 368 365 366 END IF ! dysl 369 367 CALL trace_end("compute_caldyn_vert_nh") 370 368 -
codes/icosagcm/trunk/src/dynamics/caldyn_kernels_hevi.F90
r561 r580 9 9 PRIVATE 10 10 11 REAL(rstd), PARAMETER :: pbot=1e5, Phi_bot=0., rho_bot=1e6 ! FIXME11 REAL(rstd), PARAMETER :: pbot=1e5, rho_bot=1e6 12 12 13 13 LOGICAL, SAVE :: debug_hevi_solver = .FALSE. 14 LOGICAL, PARAMETER :: rigid=.TRUE.15 14 16 15 PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_Coriolis, & … … 62 61 CALL trace_start("compute_pvort_only") 63 62 !!! Compute shallow-water potential vorticity 64 #ifdef CPP_DYSL 63 IF(dysl_pvort_only) THEN 65 64 #include "../kernels/pvort_only.k90" 66 #else 65 ELSE 66 67 67 radius_m2=radius**(-2) 68 68 DO l = ll_begin,ll_end … … 94 94 95 95 ENDDO 96 #endif 96 97 END IF ! dysl 97 98 CALL trace_end("compute_pvort_only") 98 99 99 100 END SUBROUTINE compute_pvort_only 100 101 101 SUBROUTINE compute_NH_geopot(tau, m_ik, m_il, theta, W_il, Phi_il)102 SUBROUTINE compute_NH_geopot(tau, phis, m_ik, m_il, theta, W_il, Phi_il) 102 103 REAL(rstd),INTENT(IN) :: tau ! solve Phi-tau*dPhi/dt = Phi_rhs 104 REAL(rstd),INTENT(IN) :: phis(iim*jjm) 103 105 REAL(rstd),INTENT(IN) :: m_ik(iim*jjm,llm) 104 106 REAL(rstd),INTENT(IN) :: m_il(iim*jjm,llm+1) … … 124 126 ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1 125 127 126 #ifdef CPP_DYSL 128 IF(dysl) THEN 129 #define PHI_BOT(ij) phis(ij) 127 130 #include "../kernels/compute_NH_geopot.k90" 128 #else 131 #undef PHI_BOT 132 ELSE 129 133 ! FIXME : vertical OpenMP parallelism will not work 130 134 … … 164 168 B_il(ij,1) = A_ik(ij,1) + ml_g2 + tau2_g*rho_bot 165 169 R_il(ij,1) = ml_g2*( Phi_il(ij,1)-Phi_star_il(ij,1)) & 166 + tau2_g*( p_ik(ij,1)-pbot+rho_bot*(Phi_il(ij,1)- Phi_bot) )170 + tau2_g*( p_ik(ij,1)-pbot+rho_bot*(Phi_il(ij,1)-phis(ij)) ) 167 171 ENDDO 168 172 ! inner interfaces … … 260 264 END DO ! Newton-Raphson 261 265 262 #endif 266 END IF ! dysl 263 267 264 268 END SUBROUTINE compute_NH_geopot 265 269 266 SUBROUTINE compute_caldyn_solver(tau, rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du)270 SUBROUTINE compute_caldyn_solver(tau,phis, rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du) 267 271 REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs 272 REAL(rstd),INTENT(IN) :: phis(iim*jjm) 268 273 REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 269 274 REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm,nqdyn) … … 278 283 279 284 REAL(rstd) :: berni(iim*jjm,llm) ! (W/m_il)^2 285 REAL(rstd) :: berni1(iim*jjm) ! (W/m_il)^2 280 286 REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Rd, Cvd 281 287 INTEGER :: ij, l … … 285 291 Rd=cpp*kappa 286 292 287 #ifdef CPP_DYSL 293 IF(dysl) THEN 294 288 295 !$OMP BARRIER 296 #define PHI_BOT(ij) phis(ij) 297 #define PHI_BOT_VAR phis 289 298 #include "../kernels/caldyn_solver.k90" 299 #undef PHI_BOT_VAR 300 #undef PHI_BOT 290 301 !$OMP BARRIER 291 #else 292 #define BERNI(ij) berni(ij,1) 302 303 ELSE 304 305 #define BERNI(ij) berni1(ij) 293 306 ! FIXME : vertical OpenMP parallelism will not work 294 307 … … 310 323 311 324 IF(tau>0) THEN ! solve implicit problem for geopotential 312 CALL compute_NH_geopot(tau, rhodz, m_il, theta, W, geopot)325 CALL compute_NH_geopot(tau, phis, rhodz, m_il, theta, W, geopot) 313 326 END IF 314 327 … … 372 385 ENDDO 373 386 #undef BERNI 374 #endif 387 388 END IF ! dysl 375 389 376 390 CALL trace_end("compute_caldyn_solver") … … 396 410 Rd=cpp*kappa 397 411 398 #ifdef CPP_DYSL 412 IF(dysl_caldyn_fast) THEN 399 413 #include "../kernels/caldyn_fast.k90" 400 #else 414 ELSE 415 401 416 ! Compute Bernoulli term 402 417 IF(boussinesq) THEN … … 496 511 END IF 497 512 END DO 498 #endif 513 514 END IF ! dysl 499 515 CALL trace_end("compute_caldyn_fast") 500 516 … … 515 531 CALL trace_start("compute_caldyn_Coriolis") 516 532 517 #ifdef CPP_DYSL 533 IF(dysl_caldyn_coriolis) THEN 534 518 535 #include "../kernels/coriolis.k90" 519 #else 536 537 ELSE 520 538 #define FTHETA(ij) Ftheta(ij,1) 521 539 … … 653 671 END SELECT 654 672 #undef FTHETA 655 #endif 673 674 END IF ! dysl 656 675 657 676 CALL trace_end("compute_caldyn_Coriolis") … … 667 686 668 687 REAL(rstd) :: berni(iim*jjm,llm) ! Bernoulli function 688 REAL(rstd) :: berni1(iim*jjm) ! Bernoulli function 669 689 REAL(rstd) :: uu_right, uu_lup, uu_ldown, ke, uu 670 690 INTEGER :: ij,l … … 672 692 CALL trace_start("compute_caldyn_slow_hydro") 673 693 674 #ifdef CPP_DYSL 694 IF(dysl_slow_hydro) THEN 695 675 696 #define BERNI(ij,l) berni(ij,l) 676 697 #include "../kernels/caldyn_slow_hydro.k90" 677 698 #undef BERNI 678 #else 679 #define BERNI(ij) berni(ij,1) 699 700 ELSE 701 702 #define BERNI(ij) berni1(ij) 680 703 681 704 DO l = ll_begin, ll_end … … 725 748 END IF 726 749 END DO 750 727 751 #undef BERNI 728 #endif 752 753 END IF ! dysl 729 754 CALL trace_end("compute_caldyn_slow_hydro") 730 755 END SUBROUTINE compute_caldyn_slow_hydro … … 749 774 REAL(rstd) :: W_el, W2_el, uu_right, uu_lup, uu_ldown, gPhi2, dP, divG, u2, uu 750 775 751 #ifdef CPP_DYSL752 776 REAL(rstd) :: berni(iim*jjm,llm) ! Bernoulli function 753 777 REAL(rstd) :: G_el(3*iim*jjm,llm+1) ! horizontal flux of W 754 778 REAL(rstd) :: v_el(3*iim*jjm,llm+1) 755 #else 756 REAL(rstd) :: berni(iim*jjm) ! Bernoulli function 757 REAL(rstd) :: G_el(3*iim*jjm) ! horizontal flux of W 758 REAL(rstd) :: v_el(3*iim*jjm) 759 #endif 779 780 REAL(rstd) :: berni1(iim*jjm) ! Bernoulli function 781 REAL(rstd) :: G_el1(3*iim*jjm) ! horizontal flux of W 782 REAL(rstd) :: v_el1(3*iim*jjm) 760 783 761 784 CALL trace_start("compute_caldyn_slow_NH") 762 785 763 #ifdef CPP_DYSL 786 IF(dysl) THEN 787 764 788 !$OMP BARRIER 765 789 #include "../kernels/caldyn_slow_NH.k90" 766 790 !$OMP BARRIER 767 #else 791 792 ELSE 793 794 #define BERNI(ij) berni1(ij) 795 #define G_EL(ij) G_el1(ij) 796 #define V_EL(ij) v_el1(ij) 797 768 798 DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces) 769 799 IF(l==1) THEN … … 792 822 W2_el = .5*le_de(ij+u_right) * & 793 823 ( W(ij,l)*w_il(ij,l) + W(ij+t_right,l)*w_il(ij+t_right,l) ) 794 v_el(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked795 G_ el(ij+u_right) = v_el(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el824 V_EL(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked 825 G_EL(ij+u_right) = V_EL(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el 796 826 ! Compute on edge 'lup' 797 827 W_el = .5*( W(ij,l)+W(ij+t_lup,l) ) … … 800 830 W2_el = .5*le_de(ij+u_lup) * & 801 831 ( W(ij,l)*w_il(ij,l) + W(ij+t_lup,l)*w_il(ij+t_lup,l) ) 802 v_el(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked803 G_ el(ij+u_lup) = v_el(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el832 V_EL(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked 833 G_EL(ij+u_lup) = V_EL(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el 804 834 ! Compute on edge 'ldown' 805 835 W_el = .5*( W(ij,l)+W(ij+t_ldown,l) ) … … 808 838 W2_el = .5*le_de(ij+u_ldown) * & 809 839 ( W(ij,l)*w_il(ij,l) + W(ij+t_ldown,l)*w_il(ij+t_ldown,l) ) 810 v_el(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked811 G_ el(ij+u_ldown) = v_el(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el840 V_EL(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked 841 G_EL(ij+u_ldown) = V_EL(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el 812 842 END DO 813 843 ! compute GradPhi2, dPhi, dW … … 820 850 le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)**2 + & 821 851 le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)**2 ) 822 ! gradPhi2(ij,l) = 0. ! FIXME !!823 852 824 853 dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l) -1/(2*Ai(ij))* & 825 ( DePhil(ij+u_right,l)* v_el(ij+u_right) + & ! -v.gradPhi,826 DePhil(ij+u_rup,l)* v_el(ij+u_rup) + & ! v_el already has le_de827 DePhil(ij+u_lup,l)* v_el(ij+u_lup) + &828 DePhil(ij+u_left,l)* v_el(ij+u_left) + &829 DePhil(ij+u_ldown,l)* v_el(ij+u_ldown) + &830 DePhil(ij+u_rdown,l)* v_el(ij+u_rdown) )854 ( DePhil(ij+u_right,l)*V_EL(ij+u_right) + & ! -v.gradPhi, 855 DePhil(ij+u_rup,l)*V_EL(ij+u_rup) + & ! v_el already has le_de 856 DePhil(ij+u_lup,l)*V_EL(ij+u_lup) + & 857 DePhil(ij+u_left,l)*V_EL(ij+u_left) + & 858 DePhil(ij+u_ldown,l)*V_EL(ij+u_ldown) + & 859 DePhil(ij+u_rdown,l)*V_EL(ij+u_rdown) ) 831 860 832 861 dW(ij,l) = -1./Ai(ij)*( & ! -div(G_el), 833 ne_right*G_ el(ij+u_right) + & ! G_el already has le_de834 ne_rup*G_ el(ij+u_rup) + &835 ne_lup*G_ el(ij+u_lup) + &836 ne_left*G_ el(ij+u_left) + &837 ne_ldown*G_ el(ij+u_ldown) + &838 ne_rdown*G_ el(ij+u_rdown))862 ne_right*G_EL(ij+u_right) + & ! G_el already has le_de 863 ne_rup*G_EL(ij+u_rup) + & 864 ne_lup*G_EL(ij+u_lup) + & 865 ne_left*G_EL(ij+u_left) + & 866 ne_ldown*G_EL(ij+u_ldown) + & 867 ne_rdown*G_EL(ij+u_rdown)) 839 868 END DO 840 869 END DO … … 843 872 ! Compute berni at scalar points 844 873 DO ij=ij_begin_ext, ij_end_ext 845 berni(ij) = &874 BERNI(ij) = & 846 875 1/(4*Ai(ij))*( & 847 876 le_de(ij+u_right)*u(ij+u_right,l)**2 + & … … 860 889 -0.5*(F_el(ij+u_right,l)+F_el(ij+u_right,l+1)) 861 890 hflux(ij+u_right,l) = uu_right*le_de(ij+u_right) 862 du(ij+u_right,l) = ne_right*( berni(ij)-berni(ij+t_right))891 du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right)) 863 892 ! Compute on edge 'lup' 864 893 uu_lup = 0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) & 865 894 -0.5*(F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1)) 866 895 hflux(ij+u_lup,l) = uu_lup*le_de(ij+u_lup) 867 du(ij+u_lup,l) = ne_lup*( berni(ij)-berni(ij+t_lup))896 du(ij+u_lup,l) = ne_lup*(BERNI(ij)-BERNI(ij+t_lup)) 868 897 ! Compute on edge 'ldown' 869 898 uu_ldown = 0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) & 870 899 -0.5*(F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1)) 871 900 hflux(ij+u_ldown,l) = uu_ldown*le_de(ij+u_ldown) 872 du(ij+u_ldown,l) = ne_ldown*( berni(ij)-berni(ij+t_ldown))901 du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown)) 873 902 END DO 874 903 END DO 875 #endif 904 905 #undef V_EL 906 #undef G_EL 907 #undef BERNI 908 909 END IF ! dysl 876 910 877 911 CALL trace_end("compute_caldyn_slow_NH")
Note: See TracChangeset
for help on using the changeset viewer.