![(please configure the [header_logo] section in trac.ini)](/dynamico/chrome/site/your_project_logo.png){kind=logo}
Changeset 529 for codes/icosagcm/devel/src/caldyn_kernels.f90
- Timestamp:
- 01/27/17 16:54:36 (7 years ago)
- File:
-
- 1 edited
-
codes/icosagcm/devel/src/caldyn_kernels.f90 (modified) (19 diffs)
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/caldyn_kernels.f90
r428 r529 32 32 SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv) 33 33 USE icosa 34 USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass 34 USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass, ptop 35 35 USE trace 36 36 USE omp_para … … 48 48 CALL trace_start("compute_pvort") 49 49 50 IF(caldyn_eta==eta_mass) THEN51 CALL wait_message(req_ps) ! COM0052 ELSE53 CALL wait_message(req_mass) ! COM0054 END IF55 CALL wait_message(req_theta_rhodz) ! COM0150 ! IF(caldyn_eta==eta_mass) THEN 51 ! CALL wait_message(req_ps) ! COM00 52 ! ELSE 53 ! CALL wait_message(req_mass) ! COM00 54 ! END IF 55 ! CALL wait_message(req_theta_rhodz) ! COM01 56 56 57 57 IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta 58 58 DO l = ll_begin,ll_end 59 CALL test_message(req_u)59 ! CALL test_message(req_u) 60 60 !DIR$ SIMD 61 61 DO ij=ij_begin_ext,ij_end_ext 62 m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g63 rhodz(ij,l) = m 62 m = mass_dak(l)+(ps(ij)*g+ptop)*mass_dbk(l) ! ps is actually Ms 63 rhodz(ij,l) = m/g 64 64 theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 65 65 ENDDO … … 67 67 ELSE ! Compute only theta 68 68 DO l = ll_begin,ll_end 69 CALL test_message(req_u)69 ! CALL test_message(req_u) 70 70 !DIR$ SIMD 71 71 DO ij=ij_begin_ext,ij_end_ext … … 75 75 END IF 76 76 77 CALL wait_message(req_u) ! COM0277 ! CALL wait_message(req_u) ! COM02 78 78 79 79 !!! Compute shallow-water potential vorticity … … 81 81 !DIR$ SIMD 82 82 DO ij=ij_begin_ext,ij_end_ext 83 etav= 1./Av(ij+z_up)*( ne_rup * u(ij+u_rup,l) * de(ij+u_rup) & 84 + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left) & 85 - ne_lup * u(ij+u_lup,l) * de(ij+u_lup) ) 86 87 hv = Riv2(ij,vup) * rhodz(ij,l) & 88 + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l) & 83 etav= 1./Av(ij+z_up)*( ne_rup * u(ij+u_rup,l) & 84 + ne_left * u(ij+t_rup+u_left,l) & 85 - ne_lup * u(ij+u_lup,l) ) 86 hv = Riv2(ij,vup) * rhodz(ij,l) & 87 + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l) & 89 88 + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l) 90 91 89 qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv 92 93 etav = 1./Av(ij+z_down)*( ne_ldown * u(ij+u_ldown,l) * de(ij+u_ldown) & 94 + ne_right * u(ij+t_ldown+u_right,l) * de(ij+t_ldown+u_right) & 95 - ne_rdown * u(ij+u_rdown,l) * de(ij+u_rdown) ) 96 97 hv = Riv2(ij,vdown) * rhodz(ij,l) & 90 91 etav = 1./Av(ij+z_down)*( ne_ldown * u(ij+u_ldown,l) & 92 + ne_right * u(ij+t_ldown+u_right,l) & 93 - ne_rdown * u(ij+u_rdown,l) ) 94 hv = Riv2(ij,vdown) * rhodz(ij,l) & 98 95 + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l) & 99 96 + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l) 100 101 97 qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv 102 103 98 ENDDO 104 99 … … 115 110 END SUBROUTINE compute_pvort 116 111 117 !************** *********** caldyn_horiz = caldyn_fast + caldyn_slow **********************112 !************** caldyn_horiz = caldyn_fast + caldyn_slow + caldyn_coriolis *************** 118 113 119 114 SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,convm, dtheta_rhodz, du) … … 139 134 REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux 140 135 REAL(rstd) :: berni(iim*jjm,llm) ! Bernoulli function 136 REAL(rstd) :: uu_right, uu_lup, uu_ldown 141 137 142 138 INTEGER :: i,j,ij,l … … 149 145 DO l = ll_begin, ll_end 150 146 !!! Compute mass and theta fluxes 151 IF (caldyn_conserv==energy) CALL test_message(req_qu)147 ! IF (caldyn_conserv==energy) CALL test_message(req_qu) 152 148 !DIR$ SIMD 153 149 DO ij=ij_begin_ext,ij_end_ext 154 hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) 155 hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) 156 hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown) 157 150 uu_right=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) 151 uu_lup=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) 152 uu_ldown=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) 153 uu_right= uu_right*le_de(ij+u_right) 154 uu_lup = uu_lup *le_de(ij+u_lup) 155 uu_ldown= uu_ldown*le_de(ij+u_ldown) 156 hflux(ij+u_right,l)=uu_right 157 hflux(ij+u_lup,l) =uu_lup 158 hflux(ij+u_ldown,l)=uu_ldown 159 ! hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) 160 ! hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) 161 ! hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown) 158 162 Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*hflux(ij+u_right,l) 159 163 Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l) … … 189 193 CASE(energy) ! energy-conserving TRiSK 190 194 191 CALL wait_message(req_qu) ! COM03195 ! CALL wait_message(req_qu) ! COM03 192 196 193 197 DO l=ll_begin,ll_end … … 205 209 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))+ & 206 210 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)) 207 du(ij+u_right,l) = .5*uu /de(ij+u_right)211 du(ij+u_right,l) = .5*uu 208 212 209 213 uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) + & … … 217 221 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)) + & 218 222 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)) 219 du(ij+u_lup,l) = .5*uu/de(ij+u_lup) 220 223 du(ij+u_lup,l) = .5*uu 221 224 222 225 uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) + & … … 230 233 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)) + & 231 234 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)) 232 du(ij+u_ldown,l) = .5*uu /de(ij+u_ldown)235 du(ij+u_ldown,l) = .5*uu 233 236 234 237 ENDDO … … 251 254 wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+ & 252 255 wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l) 253 du(ij+u_right,l) = qu(ij+u_right,l)*uu /de(ij+u_right)256 du(ij+u_right,l) = qu(ij+u_right,l)*uu 254 257 255 258 … … 264 267 wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+ & 265 268 wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l) 266 du(ij+u_lup,l) = qu(ij+u_lup,l)*uu /de(ij+u_lup)269 du(ij+u_lup,l) = qu(ij+u_lup,l)*uu 267 270 268 271 uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+ & … … 276 279 wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+ & 277 280 wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l) 278 du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu /de(ij+u_ldown)281 du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu 279 282 280 283 ENDDO … … 290 293 !DIR$ SIMD 291 294 DO ij=ij_begin,ij_end 292 293 berni(ij,l) = pk(ij,l) + & 294 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 + & 295 le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 + & 296 le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 + & 297 le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 + & 298 le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 + & 299 le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 295 berni(ij,l) = pk(ij,l) + 1/(4*Ai(ij))*( & 296 le_de(ij+u_right)*u(ij+u_right,l)**2 + & 297 le_de(ij+u_rup)*u(ij+u_rup,l)**2 + & 298 le_de(ij+u_lup)*u(ij+u_lup,l)**2 + & 299 le_de(ij+u_left)*u(ij+u_left,l)**2 + & 300 le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 + & 301 le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 300 302 ! from now on pk contains the vertically-averaged geopotential 301 303 pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) … … 308 310 !DIR$ SIMD 309 311 DO ij=ij_begin,ij_end 310 311 312 berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 312 + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 + & 313 le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 + & 314 le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 + & 315 le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 + & 316 le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 + & 317 le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 313 + 1/(4*Ai(ij))*( & 314 le_de(ij+u_right)*u(ij+u_right,l)**2 + & 315 le_de(ij+u_rup)*u(ij+u_rup,l)**2 + & 316 le_de(ij+u_lup)*u(ij+u_lup,l)**2 + & 317 le_de(ij+u_left)*u(ij+u_left,l)**2 + & 318 le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 + & 319 le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 318 320 ENDDO 319 321 ENDDO … … 326 328 DO ij=ij_begin,ij_end 327 329 328 du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) *( &330 du(ij+u_right,l) = du(ij+u_right,l) + ( & 329 331 0.5*(theta(ij,l)+theta(ij+t_right,l)) & 330 332 *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l)) & … … 332 334 333 335 334 du(ij+u_lup,l) = du(ij+u_lup,l) + 1/de(ij+u_lup) *( &336 du(ij+u_lup,l) = du(ij+u_lup,l) + ( & 335 337 0.5*(theta(ij,l)+theta(ij+t_lup,l)) & 336 338 *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l)) & 337 339 + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) ) 338 340 339 du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) *( &341 du(ij+u_ldown,l) = du(ij+u_ldown,l) + ( & 340 342 0.5*(theta(ij,l)+theta(ij+t_ldown,l)) & 341 343 *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l)) &
Note: See TracChangeset
for help on using the changeset viewer.
