Changeset 734 for codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90
- Timestamp:
- 08/27/18 13:51:26 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90
r733 r734 13 13 LOGICAL, SAVE :: debug_hevi_solver = .FALSE. 14 14 15 PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_ Coriolis, &15 PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_Kv, compute_caldyn_Coriolis, & 16 16 compute_caldyn_slow_hydro, compute_caldyn_slow_NH, & 17 17 compute_caldyn_solver, compute_caldyn_fast … … 99 99 100 100 END SUBROUTINE compute_pvort_only 101 102 SUBROUTINE compute_caldyn_kv(ue, Kv) 103 REAL(rstd),INTENT(IN) :: ue(3*iim*jjm,llm) 104 REAL(rstd),INTENT(OUT) :: Kv(2*iim*jjm,llm) 105 REAL(rstd) :: ue2(3*iim*jjm), dem2(3*iim*jjm), r2_Av(2*iim*jjm), rad2 106 INTEGER :: ij,l, u_up, u_down 107 108 u_up = t_lup + u_right 109 u_down = t_rdown + u_left 110 111 rad2=radius**2 112 113 !DIR$ SIMD 114 DO ij=ij_begin_ext,ij_end_ext 115 dem2(ij+u_right) = de(ij+u_right)**(-2) 116 dem2(ij+u_lup) = de(ij+u_lup)**(-2) 117 dem2(ij+u_ldown) = de(ij+u_ldown)**(-2) 118 r2_Av(ij+z_up) = rad2*(1./Av(ij+z_up)) 119 r2_Av(ij+z_down) = rad2*(1./Av(ij+z_down)) 120 END DO 121 122 DO l=ll_begin,ll_end 123 ! compute squared normal component from 1-form 124 !DIR$ SIMD 125 DO ij=ij_begin_ext,ij_end_ext 126 ue2(ij+u_right) = dem2(ij+u_right)* (ue(ij+u_right,l)**2) 127 ue2(ij+u_lup) = dem2(ij+u_lup) * (ue(ij+u_lup,l)**2) 128 ue2(ij+u_ldown) = dem2(ij+u_ldown)* (ue(ij+u_ldown,l)**2) 129 END DO 130 ! average squared normal component to vertices 131 !DIR$ SIMD 132 DO ij=ij_begin_ext,ij_end_ext 133 Kv(ij+z_up,l) = r2_Av(ij+z_up)*( & 134 S1(ij,vup)*ue2(ij+u_rup) + & 135 S2(ij,vup)*ue2(ij+u_lup) + & 136 S2(ij+t_lup,vrdown)*ue2(ij+u_up)) 137 138 Kv(ij+z_down,l) = r2_Av(ij+z_down)*( & 139 S1(ij,vdown)*ue2(ij+u_ldown) + & 140 S2(ij,vdown)*ue2(ij+u_rdown) + & 141 S2(ij+t_rdown,vlup)*ue2(ij+u_down) ) 142 ENDDO 143 ENDDO 144 END SUBROUTINE compute_caldyn_kv 101 145 102 146 SUBROUTINE compute_NH_geopot(tau, phis, m_ik, m_il, theta, W_il, Phi_il) … … 663 707 wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l) 664 708 665 du(ij+u_right,l) = du(ij+u_right,l) + .5*uu_right666 du(ij+u_lup,l) = du(ij+u_lup,l) + .5*uu_lup667 du(ij+u_ldown,l) = du(ij+u_ldown,l) + .5*uu_ldown709 du(ij+u_right,l) = du(ij+u_right,l) + uu_right*qu(ij+u_right,l) 710 du(ij+u_lup,l) = du(ij+u_lup,l) + uu_lup*qu(ij+u_lup,l) 711 du(ij+u_ldown,l) = du(ij+u_ldown,l) + uu_ldown*qu(ij+u_ldown,l) 668 712 END DO 669 713 END DO … … 679 723 END SUBROUTINE compute_caldyn_Coriolis 680 724 681 SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux, du, zero)725 SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,Kv,du, zero) 682 726 LOGICAL, INTENT(IN) :: zero 683 727 REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm) ! prognostic "velocity" 728 REAL(rstd),INTENT(IN) :: Kv(2*iim*jjm,llm) ! kinetic energy at vertices 684 729 REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 685 730 REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s … … 719 764 ENDDO 720 765 ! Compute Bernoulli=kinetic energy 721 !DIR$ SIMD 722 DO ij=ij_begin,ij_end 723 BERNI(ij) = & 724 1/(4*Ai(ij))*(le_de(ij+u_right)*u(ij+u_right,l)**2 + & 725 le_de(ij+u_rup)*u(ij+u_rup,l)**2 + & 726 le_de(ij+u_lup)*u(ij+u_lup,l)**2 + & 727 le_de(ij+u_left)*u(ij+u_left,l)**2 + & 728 le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 + & 729 le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 730 ENDDO 766 IF(caldyn_kinetic==kinetic_trisk) THEN 767 !DIR$ SIMD 768 DO ij=ij_begin,ij_end 769 BERNI(ij) = & 770 1/(4*Ai(ij))*(le_de(ij+u_right)*u(ij+u_right,l)**2 + & 771 le_de(ij+u_rup)*u(ij+u_rup,l)**2 + & 772 le_de(ij+u_lup)*u(ij+u_lup,l)**2 + & 773 le_de(ij+u_left)*u(ij+u_left,l)**2 + & 774 le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 + & 775 le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 776 ENDDO 777 ELSE 778 !DIR$ SIMD 779 DO ij=ij_begin,ij_end 780 BERNI(ij) = Riv(ij,vup) *Kv(ij+z_up,l) + & 781 Riv(ij,vlup) *Kv(ij+z_lup,l) + & 782 Riv(ij,vldown)*Kv(ij+z_ldown,l) + & 783 Riv(ij,vdown) *Kv(ij+z_down,l) + & 784 Riv(ij,vrdown)*Kv(ij+z_rdown,l) + & 785 Riv(ij,vrup) *Kv(ij+z_rup,l) 786 END DO 787 END IF 731 788 ! Compute du=-grad(Bernoulli) 732 789 IF(zero) THEN
Note: See TracChangeset
for help on using the changeset viewer.