Changeset 529 for codes/icosagcm/devel/src
- Timestamp:
- 01/27/17 16:54:36 (7 years ago)
- Location:
- codes/icosagcm/devel/src
- Files:
-
- 10 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/caldyn_gcm.f90
r413 r529 2 2 USE icosa 3 3 USE transfert_mod 4 USE caldyn_kernels_hevi_mod 4 5 USE caldyn_kernels_base_mod 5 6 USE caldyn_kernels_mod … … 190 191 191 192 ! temporary shared variable 192 REAL(rstd),POINTER :: theta(:,: )193 REAL(rstd),POINTER :: theta(:,:,:) 193 194 REAL(rstd),POINTER :: pk(:,:) 194 195 REAL(rstd),POINTER :: geopot(:,:) … … 225 226 IF(caldyn_eta==eta_mass) THEN 226 227 CALL send_message(f_ps,req_ps) ! COM00 228 CALL wait_message(req_ps) ! COM00 227 229 ELSE 228 230 CALL send_message(f_mass,req_mass) ! COM00 231 CALL wait_message(req_mass) ! COM00 229 232 END IF 230 233 231 234 CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01 235 CALL wait_message(req_theta_rhodz) ! COM01 Moved from caldyn_pvort 232 236 CALL send_message(f_u,req_u) ! COM02 237 CALL wait_message(req_u) ! COM02 238 239 IF(.NOT.hydrostatic) THEN 240 STOP 'caldyn_gcm may not be used yet when non-hydrostatic' 241 END IF 233 242 234 243 SELECT CASE(caldyn_conserv) … … 245 254 qu=f_qu(ind) 246 255 qv=f_qv(ind) 256 pk = f_pk(ind) 257 geopot = f_geopot(ind) 258 hflux=f_hflux(ind) 259 convm = f_dmass(ind) 260 dtheta_rhodz=f_dtheta_rhodz(ind) 261 du=f_du(ind) 247 262 CALL compute_pvort(ps,u,theta_rhodz(:,:,1), mass,theta,qu,qv) ! COM00 COM01 COM02 248 ENDDO 249 ! CALL checksum(f_mass) 250 ! CALL checksum(f_theta) 251 252 CALL send_message(f_qu,req_qu) ! COM03 wait_message in caldyn_horiz 253 ! CALL wait_message(req_qu) 263 ! CALL compute_theta(ps,theta_rhodz, mass,theta) 264 ! CALL compute_pvort_only(u,mass,qu,qv) 265 266 CALL compute_geopot(mass,theta, ps,pk,geopot) 267 ! du(:,:)=0. 268 ! CALL compute_caldyn_fast(0.,u,mass,theta,pk,geopot,du) 269 ENDDO 270 271 CALL send_message(f_u,req_u) ! COM02 272 CALL wait_message(req_u) ! COM02 273 CALL send_message(f_qu,req_qu) ! COM03 274 CALL wait_message(req_qu) ! COM03 254 275 255 276 DO ind=1,ndomain … … 259 280 ps=f_ps(ind) 260 281 u=f_u(ind) 282 theta_rhodz = f_theta_rhodz(ind) 261 283 mass=f_mass(ind) 262 284 theta = f_theta(ind) 263 285 qu=f_qu(ind) 286 qv=f_qv(ind) 264 287 pk = f_pk(ind) 265 288 geopot = f_geopot(ind) 266 CALL compute_geopot(ps,mass,theta, pk,geopot)267 289 hflux=f_hflux(ind) 268 290 convm = f_dmass(ind) 269 291 dtheta_rhodz=f_dtheta_rhodz(ind) 270 292 du=f_du(ind) 293 271 294 CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz(:,:,1),du) 295 ! CALL compute_caldyn_slow_hydro(u,mass,hflux,du, .FALSE.) ! FIXME 296 ! CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du) 272 297 IF(caldyn_eta==eta_mass) THEN 273 298 wflux=f_wflux(ind) … … 278 303 ENDDO 279 304 280 ! CALL checksum(f_geopot)281 ! CALL checksum(f_dmass)282 ! CALL checksum(f_pk)283 ! CALL checksum(f_pk)284 285 305 CASE(enstrophy) ! enstrophy-conserving 286 306 DO ind=1,ndomain -
codes/icosagcm/devel/src/caldyn_hevi.f90
r404 r529 155 155 156 156 IF(hydrostatic) THEN 157 CALL compute_caldyn_slow_hydro(u,mass,hflux,du )157 CALL compute_caldyn_slow_hydro(u,mass,hflux,du, .TRUE.) 158 158 ELSE 159 159 W = f_W(ind) -
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)) & -
codes/icosagcm/devel/src/caldyn_kernels_base.f90
r479 r529 7 7 IMPLICIT NONE 8 8 PRIVATE 9 10 LOGICAL, PARAMETER, PUBLIC :: DEC = .TRUE.11 9 12 10 INTEGER, PARAMETER,PUBLIC :: energy=1, enstrophy=2 … … 219 217 DO ij=ij_begin,ij_end 220 218 ! dps/dt = -int(div flux)dz 221 IF(DEC) THEN 222 dps(ij) = convm(ij,1) 223 ELSE 224 dps(ij) = convm(ij,1) * g 225 END IF 219 dps(ij) = convm(ij,1) 226 220 ENDDO 227 221 ENDIF -
codes/icosagcm/devel/src/caldyn_kernels_hevi.f90
r499 r529 32 32 !DIR$ SIMD 33 33 DO ij=ij_begin_ext,ij_end_ext 34 IF(DEC) THEN ! ps is actually Ms 35 m = mass_dak(l)+(ps(ij)*g+ptop)*mass_dbk(l) 36 ELSE 37 m = mass_dak(l)+ps(ij)*mass_dbk(l) 38 END IF 34 m = mass_dak(l)+(ps(ij)*g+ptop)*mass_dbk(l) ! ps is actually Ms 39 35 rhodz(ij,l) = m/g 40 36 END DO … … 69 65 !DIR$ SIMD 70 66 DO ij=ij_begin_ext,ij_end_ext 71 IF(DEC) THEN 72 etav= 1./Av(ij+z_up)*( ne_rup * u(ij+u_rup,l) & 73 + ne_left * u(ij+t_rup+u_left,l) & 74 - ne_lup * u(ij+u_lup,l) ) 75 76 hv = Riv2(ij,vup) * rhodz(ij,l) & 77 + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l) & 78 + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l) 79 qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv 80 81 etav = 1./Av(ij+z_down)*( ne_ldown * u(ij+u_ldown,l) & 82 + ne_right * u(ij+t_ldown+u_right,l) & 83 - ne_rdown * u(ij+u_rdown,l) ) 84 hv = Riv2(ij,vdown) * rhodz(ij,l) & 85 + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l) & 86 + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l) 87 qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv 88 ELSE 89 etav= 1./Av(ij+z_up)*( ne_rup * u(ij+u_rup,l) * de(ij+u_rup) & 90 + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left) & 91 - ne_lup * u(ij+u_lup,l) * de(ij+u_lup) ) 92 93 hv = Riv2(ij,vup) * rhodz(ij,l) & 94 + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l) & 95 + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l) 96 qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv 97 98 etav = 1./Av(ij+z_down)*( ne_ldown * u(ij+u_ldown,l) * de(ij+u_ldown) & 99 + ne_right * u(ij+t_ldown+u_right,l) * de(ij+t_ldown+u_right) & 100 - ne_rdown * u(ij+u_rdown,l) * de(ij+u_rdown) ) 101 hv = Riv2(ij,vdown) * rhodz(ij,l) & 102 + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l) & 103 + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l) 104 qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv 105 END IF 67 etav= 1./Av(ij+z_up)*( ne_rup * u(ij+u_rup,l) & 68 + ne_left * u(ij+t_rup+u_left,l) & 69 - ne_lup * u(ij+u_lup,l) ) 70 hv = Riv2(ij,vup) * rhodz(ij,l) & 71 + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l) & 72 + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l) 73 qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv 74 75 etav = 1./Av(ij+z_down)*( ne_ldown * u(ij+u_ldown,l) & 76 + ne_right * u(ij+t_ldown+u_right,l) & 77 - ne_rdown * u(ij+u_rdown,l) ) 78 hv = Riv2(ij,vdown) * rhodz(ij,l) & 79 + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l) & 80 + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l) 81 qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv 106 82 ENDDO 107 83 … … 649 625 END SUBROUTINE compute_caldyn_Coriolis 650 626 651 SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du) 627 SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du, zero) 628 LOGICAL, INTENT(IN) :: zero 652 629 REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm) ! prognostic "velocity" 653 630 REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 654 631 REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s 655 REAL(rstd),INTENT( OUT) :: du(3*iim*jjm,llm)632 REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm) 656 633 657 634 REAL(rstd) :: berni(iim*jjm) ! Bernoulli function … … 690 667 ENDDO 691 668 ! Compute du=-grad(Bernoulli) 692 !DIR$ SIMD 693 DO ij=ij_begin,ij_end 669 IF(zero) THEN 670 !DIR$ SIMD 671 DO ij=ij_begin,ij_end 694 672 du(ij+u_right,l) = ne_right*(berni(ij)-berni(ij+t_right)) 695 673 du(ij+u_lup,l) = ne_lup*(berni(ij)-berni(ij+t_lup)) 696 674 du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown)) 697 END DO 675 END DO 676 ELSE 677 !DIR$ SIMD 678 DO ij=ij_begin,ij_end 679 du(ij+u_right,l) = du(ij+u_right,l) + & 680 ne_right*(berni(ij)-berni(ij+t_right)) 681 du(ij+u_lup,l) = du(ij+u_lup,l) + & 682 ne_lup*(berni(ij)-berni(ij+t_lup)) 683 du(ij+u_ldown,l) = du(ij+u_ldown,l) + & 684 ne_ldown*(berni(ij)-berni(ij+t_ldown)) 685 END DO 686 END IF 698 687 END DO 699 688 -
codes/icosagcm/devel/src/dissip_gcm.f90
r487 r529 7 7 TYPE(t_field),POINTER,SAVE :: f_due_diss2(:) 8 8 9 TYPE(t_field),POINTER,SAVE :: f_theta(:), f_phi(:), f_pk(:), f_pks(:), f_p(:)10 9 TYPE(t_field),POINTER,SAVE :: f_dtheta_diss(:) 11 10 TYPE(t_field),POINTER,SAVE :: f_dtheta_rhodz_diss(:) … … 35 34 INTEGER, SAVE :: rayleigh_friction_type, rayleigh_shear 36 35 !$OMP THREADPRIVATE(rayleigh_friction_type) 36 !$OMP THREADPRIVATE(rayleigh_shear) 37 37 REAL, SAVE :: rayleigh_tau 38 !$OMP THREADPRIVATE(rayleigh_ shear)38 !$OMP THREADPRIVATE(rayleigh_tau) 39 39 40 40 REAL,SAVE :: dtdissip … … 49 49 CALL allocate_field(f_due_diss1,field_u,type_real,llm) 50 50 CALL allocate_field(f_due_diss2,field_u,type_real,llm) 51 CALL allocate_field(f_theta,field_t,type_real,llm)52 51 CALL allocate_field(f_dtheta_diss,field_t,type_real,llm) 53 52 CALL allocate_field(f_dtheta_rhodz_diss,field_t,type_real,llm) 54 55 CALL allocate_field(f_phi,field_t,type_real,llm)56 CALL allocate_field(f_pk,field_t,type_real,llm)57 CALL allocate_field(f_p,field_t,type_real,llm+1)58 CALL allocate_field(f_pks,field_t,type_real)59 60 53 ALLOCATE(tau_graddiv(llm)) 61 54 ALLOCATE(tau_gradrot(llm)) … … 435 428 mintau=MIN(mintau,tau_divgrad(l)) 436 429 ENDDO 437 430 431 IF(rayleigh_friction_type>0) mintau=MIN(mintau,rayleigh_tau) 432 438 433 IF(mintau>0) THEN 439 434 itau_dissip=INT(mintau/dt) … … 444 439 END IF 445 440 itau_dissip=MAX(1,itau_dissip) 446 IF (is_master) PRINT *,"mintau ",mintau,"itau_dissip",itau_dissip," dtdissip ",dtdissip 441 IF (is_master) PRINT *,"rayleigh_tau",rayleigh_tau, "mintau ",mintau, & 442 "itau_dissip",itau_dissip," dtdissip ",dtdissip 447 443 448 444 END SUBROUTINE init_dissip 449 445 450 446 451 SUBROUTINE dissip(f_ ue,f_due,f_mass,f_phis,f_theta_rhodz,f_dtheta_rhodz)447 SUBROUTINE dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_ue, f_dtheta_rhodz,f_due) 452 448 USE icosa 453 449 USE theta2theta_rhodz_mod … … 459 455 USE omp_para 460 456 IMPLICIT NONE 461 TYPE(t_field),POINTER :: f_ue(:) 462 TYPE(t_field),POINTER :: f_due(:) 463 TYPE(t_field),POINTER :: f_mass(:), f_phis(:) 464 TYPE(t_field),POINTER :: f_theta_rhodz(:) 465 TYPE(t_field),POINTER :: f_dtheta_rhodz(:) 457 TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_phis(:), f_geopot(:) 458 TYPE(t_field),POINTER :: f_theta_rhodz(:), f_dtheta_rhodz(:) 459 TYPE(t_field),POINTER :: f_ue(:), f_due(:) 466 460 467 461 REAL(rstd),POINTER :: due(:,:) … … 483 477 CALL divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz_diss) 484 478 485 ! later for openmp486 ! IF(rayleigh_friction_type>0) THEN487 ! CALL pression(f_ps, f_p)488 ! CALL exner(f_ps, f_p, f_pks, f_pk)489 ! CALL geopotential(f_phis,f_pks,f_pk,f_theta,f_phi)490 ! END IF491 492 479 DO ind=1,ndomain 493 480 IF (.NOT. assigned_domain(ind)) CYCLE … … 515 502 ! due=0 516 503 517 ! later for openmp 518 ! IF(rayleigh_friction_type>0) THEN 519 ! phi=f_phi(ind) 520 ! ue=f_ue(ind) 521 ! DO l=1,llm 522 ! DO j=jj_begin,jj_end 523 ! DO i=ii_begin,ii_end 524 ! n=(j-1)*iim+i 525 ! CALL relax(t_right, u_right) 526 ! CALL relax(t_lup, u_lup) 527 ! CALL relax(t_ldown, u_ldown) 528 ! ENDDO 529 ! ENDDO 530 ! END DO 531 ! END IF 504 IF(rayleigh_friction_type>0) THEN 505 phi=f_geopot(ind) 506 ue=f_ue(ind) 507 DO l=ll_begin,ll_end 508 DO ij=ij_begin,ij_end 509 CALL relax(t_right, u_right) 510 CALL relax(t_lup, u_lup) 511 CALL relax(t_ldown, u_ldown) 512 ENDDO 513 END DO 514 END IF 532 515 END DO 533 516 … … 558 541 fz = sin((pi/2)*(z-zh)/(ztop-zh)) 559 542 fz = fz*fz/rayleigh_tau 560 ! fz = 1/rayleigh_tau 561 due(nn,l) = due(nn,l) - fz*(ue(nn,l)-uref) 543 due(nn,l) = due(nn,l) - itau_dissip*fz*(ue(nn,l)-uref) 544 ! fz = 1./rayleigh_tau 562 545 ! due(nn,l) = due(nn,l) - fz*ue(nn,l) 563 546 END IF -
codes/icosagcm/devel/src/explicit_scheme.f90
r487 r529 26 26 REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:), du(:,:) 27 27 REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), massm1(:,:), massm2(:,:), dmass(:,:) 28 REAL(rstd),POINTER :: theta_rhodz(:,: ) , theta_rhodzm1(:,:), theta_rhodzm2(:,:), dtheta_rhodz(:,:)28 REAL(rstd),POINTER :: theta_rhodz(:,:,:) , theta_rhodzm1(:,:,:), theta_rhodzm2(:,:,:), dtheta_rhodz(:,:,:) 29 29 REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:) 30 30 LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time 31 31 INTEGER :: it,stage 32 33 CALL legacy_to_DEC(f_ps, f_u) 32 34 DO stage=1,nb_stage 33 35 ! CALL checksum(f_ps) … … 53 55 END SELECT 54 56 END DO 57 CALL DEC_to_legacy(f_ps, f_u) 55 58 56 59 CONTAINS … … 142 145 um1(ij+u_lup,l)=u(ij+u_lup,l) 143 146 um1(ij+u_ldown,l)=u(ij+u_ldown,l) 144 theta_rhodzm1(ij,l )=theta_rhodz(ij,l)147 theta_rhodzm1(ij,l,1)=theta_rhodz(ij,l,1) 145 148 ENDDO 146 149 ENDDO … … 153 156 u(ij+u_lup,l)=um1(ij+u_lup,l)+tau*du(ij+u_lup,l) 154 157 u(ij+u_ldown,l)=um1(ij+u_ldown,l)+tau*du(ij+u_ldown,l) 155 theta_rhodz(ij,l )=theta_rhodzm1(ij,l)+tau*dtheta_rhodz(ij,l)158 theta_rhodz(ij,l,1)=theta_rhodzm1(ij,l,1)+tau*dtheta_rhodz(ij,l,1) 156 159 ENDDO 157 160 ENDDO … … 188 191 189 192 IF (stage==1) THEN 190 psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,: )=theta_rhodz(:,:)193 psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:,:)=theta_rhodz(:,:,:) 191 194 ps(:)=psm1(:)+tau*dps(:) 192 195 u(:,:)=um1(:,:)+tau*du(:,:) 193 theta_rhodz(:,: )=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)196 theta_rhodz(:,:,:)=theta_rhodzm1(:,:,:)+tau*dtheta_rhodz(:,:,:) 194 197 195 198 ELSE IF (stage==2) THEN … … 197 200 ps(:)=psm1(:)+tau*dps(:) 198 201 u(:,:)=um1(:,:)+tau*du(:,:) 199 theta_rhodz(:,: )=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)200 201 psm2(:)=psm1(:) ; theta_rhodzm2(:,: )=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)202 psm1(:)=ps(:) ; theta_rhodzm1(:,: )=theta_rhodz(:,:) ; um1(:,:)=u(:,:)202 theta_rhodz(:,:,:)=theta_rhodzm1(:,:,:)+tau*dtheta_rhodz(:,:,:) 203 204 psm2(:)=psm1(:) ; theta_rhodzm2(:,:,:)=theta_rhodzm1(:,:,:) ; um2(:,:)=um1(:,:) 205 psm1(:)=ps(:) ; theta_rhodzm1(:,:,:)=theta_rhodz(:,:,:) ; um1(:,:)=u(:,:) 203 206 204 207 ELSE … … 206 209 ps(:)=psm2(:)+2*tau*dps(:) 207 210 u(:,:)=um2(:,:)+2*tau*du(:,:) 208 theta_rhodz(:,: )=theta_rhodzm2(:,:)+2*tau*dtheta_rhodz(:,:)209 210 psm2(:)=psm1(:) ; theta_rhodzm2(:,: )=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)211 psm1(:)=ps(:) ; theta_rhodzm1(:,: )=theta_rhodz(:,:) ; um1(:,:)=u(:,:)211 theta_rhodz(:,:,:)=theta_rhodzm2(:,:,:)+2*tau*dtheta_rhodz(:,:,:) 212 213 psm2(:)=psm1(:) ; theta_rhodzm2(:,:,:)=theta_rhodzm1(:,:,:) ; um2(:,:)=um1(:,:) 214 psm1(:)=ps(:) ; theta_rhodzm1(:,:,:)=theta_rhodz(:,:,:) ; um1(:,:)=u(:,:) 212 215 213 216 ENDIF -
codes/icosagcm/devel/src/hevi_scheme.f90
r480 r529 4 4 USE field_mod 5 5 USE euler_scheme_mod 6 USE caldyn_kernels_base_mod, ONLY : DEC7 6 IMPLICIT NONE 8 7 PRIVATE … … 60 59 REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:) 61 60 62 IF(DEC)CALL legacy_to_DEC(f_ps, f_u)61 CALL legacy_to_DEC(f_ps, f_u) 63 62 DO j=1,nb_stage 64 63 CALL caldyn_hevi((j==1) .AND. (MOD(it,itau_out)==0), taujj(j), & … … 97 96 END DO 98 97 END DO 99 IF(DEC)CALL DEC_to_legacy(f_ps, f_u)98 CALL DEC_to_legacy(f_ps, f_u) 100 99 END SUBROUTINE HEVI_scheme 101 100 -
codes/icosagcm/devel/src/observable.f90
r482 r529 77 77 END IF 78 78 79 CALL divide_by_mass(1, f_mass, f_theta_rhodz, f_buf_i)80 IF(init) THEN81 CALL output_field("theta_init",f_buf_i)82 ELSE83 CALL output_field("theta",f_buf_i)84 END IF85 86 79 IF(nqdyn>1) THEN 87 80 CALL divide_by_mass(2, f_mass, f_theta_rhodz, f_buf_i) … … 93 86 END IF 94 87 95 ! CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) 96 ! CALL Tv2T(f_buf_i,f_q,f_buf1_i) 97 CALL diagnose_temperature(f_ps, f_theta_rhodz, f_q, f_buf_i) 88 CALL divide_by_mass(1, f_mass, f_theta_rhodz, f_buf_i) 89 IF(init) THEN 90 CALL output_field("theta_init",f_buf_i) 91 ELSE 92 CALL output_field("theta",f_buf_i) 93 END IF 94 95 CALL pression_mid(f_ps, f_pmid) 96 CALL diagnose_temperature(f_pmid, f_q, f_buf_i) ! f_buf_i : IN = theta, out = T 98 97 99 98 IF(init) THEN … … 112 111 CALL transfert_request(f_buf_uh,req_e1_vect) 113 112 CALL un2ulonlat(f_buf_uh, f_buf_ulon, f_buf_ulat) 114 CALL pression_mid(f_ps, f_pmid)115 113 IF(init) THEN 116 114 CALL output_field("uz_init",f_buf_i) … … 236 234 END SUBROUTINE compute_prognostic_vel_to_horiz 237 235 238 SUBROUTINE diagnose_temperature(f_p s,f_theta_rhodz,f_q,f_temp)236 SUBROUTINE diagnose_temperature(f_pmid,f_q,f_temp) 239 237 USE icosa 240 238 USE pression_mod 241 239 IMPLICIT NONE 242 TYPE(t_field), POINTER :: f_ps(:) ! IN 243 TYPE(t_field), POINTER :: f_theta_rhodz(:) ! IN 240 TYPE(t_field), POINTER :: f_pmid(:) ! IN 244 241 TYPE(t_field), POINTER :: f_q(:) ! IN 245 TYPE(t_field), POINTER :: f_temp(:) ! OUT 246 247 REAL(rstd), POINTER :: ps(:) 248 REAL(rstd), POINTER :: theta_rhodz(:,:,:) 242 TYPE(t_field), POINTER :: f_temp(:) ! INOUT 243 244 REAL(rstd), POINTER :: pmid(:,:) 249 245 REAL(rstd), POINTER :: q(:,:,:) 250 246 REAL(rstd), POINTER :: temp(:,:) … … 255 251 CALL swap_dimensions(ind) 256 252 CALL swap_geometry(ind) 257 ps=f_ps(ind) 258 theta_rhodz=f_theta_rhodz(ind) 253 pmid=f_pmid(ind) 259 254 q=f_q(ind) 260 255 temp=f_temp(ind) 261 CALL compute_diagnose_temp(p s,theta_rhodz,q,temp)256 CALL compute_diagnose_temp(pmid,q,temp) 262 257 END DO 263 258 END SUBROUTINE diagnose_temperature 264 259 265 SUBROUTINE compute_diagnose_temp(p s,theta_rhodz,q,temp)260 SUBROUTINE compute_diagnose_temp(pmid,q,temp) 266 261 USE omp_para 267 262 USE pression_mod 268 REAL(rstd),INTENT(IN) :: ps(iim*jjm) 269 REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm,nqdyn) 270 REAL(rstd),INTENT(IN) :: q(iim*jjm,llm,nqtot) 271 REAL(rstd),INTENT(OUT) :: temp(iim*jjm,llm) 272 273 REAL(rstd) :: p(iim*jjm,llm+1) 263 REAL(rstd),INTENT(IN) :: pmid(iim*jjm,llm) 264 REAL(rstd),INTENT(IN) :: q(iim*jjm,llm,nqtot) 265 REAL(rstd),INTENT(INOUT) :: temp(iim*jjm,llm) 266 274 267 REAL(rstd) :: Rd, p_ik, theta_ik, temp_ik, qv, chi, Rmix 275 268 INTEGER :: ij,l 276 269 277 270 Rd = kappa*cpp 278 CALL compute_pression(ps,p,0)279 271 DO l=ll_begin,ll_end 280 272 DO ij=ij_begin,ij_end 281 p_ik = .5*(p(ij,l)+p(ij,l+1))282 theta_ik = g*theta_rhodz(ij,l,1)/(p(ij,l)-p(ij,l+1))273 p_ik = pmid(ij,l) 274 theta_ik = temp(ij,l) 283 275 qv = q(ij,l,1) ! water vaper mixing ratio = mv/md 284 276 SELECT CASE(caldyn_thermo) -
codes/icosagcm/devel/src/timeloop_gcm.f90
r528 r529 211 211 IF(positive_theta) CALL copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q) 212 212 213 CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,itau0) 214 215 CALL trace_on 216 217 IF (xios_output) THEN ! we must call update_calendar before any XIOS output 218 IF (is_omp_master) CALL xios_update_calendar(1) 219 END IF 220 221 IF (write_start) CALL write_etat0(itau0,f_ps, f_phis,f_theta_rhodz,f_u,f_q) ! FIXME : write_start undefined 222 223 CALL write_output_fields_basic(.TRUE., f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q) 224 213 225 !$OMP MASTER 214 226 CALL SYSTEM_CLOCK(start_clock) 215 227 CALL SYSTEM_CLOCK(count_rate=rate_clock) 216 228 !$OMP END MASTER 217 218 CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,itau0)219 220 CALL trace_on221 222 IF (xios_output) THEN ! we must call update_calendar before any XIOS output223 IF (is_omp_master) CALL xios_update_calendar(1)224 END IF225 226 IF (write_start) CALL write_etat0(itau0,f_ps, f_phis,f_theta_rhodz,f_u,f_q) ! FIXME : write_start undefined227 228 CALL write_output_fields_basic(.TRUE., f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)229 229 230 230 DO it=itau0+1,itau0+itaumax … … 284 284 f_ps,f_dps,f_u,f_theta_rhodz,f_phis) 285 285 286 CALL dissip(f_ u,f_du,f_mass,f_phis, f_theta_rhodz,f_dtheta_rhodz)286 CALL dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_u, f_dtheta_rhodz,f_du) 287 287 288 288 CALL euler_scheme(.FALSE.) ! update only u, theta
Note: See TracChangeset
for help on using the changeset viewer.