![(please configure the [header_logo] section in trac.ini)](/dynamico/chrome/site/your_project_logo.png){kind=logo}
Changeset 128
- Timestamp:
- 02/04/13 11:04:04 (11 years ago)
- File:
-
- 1 edited
-
codes/icosagcm/trunk/src/caldyn_gcm.f90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/caldyn_gcm.f90
r127 r128 31 31 STOP 32 32 END SELECT 33 PRINT *, 'caldyn_conserv=',def 33 34 34 35 def='direct' … … 108 109 SELECT CASE(caldyn_conserv) 109 110 CASE(1) ! energy-conserving 110 CASE(2) ! enstrophy-conserving111 111 DO ind=1,ndomain 112 112 CALL swap_dimensions(ind) 113 113 CALL swap_geometry(ind) 114 114 ps=f_ps(ind) 115 rhodz=f_rhodz(ind) 116 p=f_p(ind) 117 qu=f_qu(ind) 118 u=f_u(ind) 119 !$OMP PARALLEL DEFAULT(SHARED) 120 CALL compute_pvort(ps, u, p,rhodz,qu) 121 !$OMP END PARALLEL 122 ENDDO 123 124 CALL transfert_request(f_qu,req_e1) 125 126 DO ind=1,ndomain 127 CALL swap_dimensions(ind) 128 CALL swap_geometry(ind) 115 129 phis=f_phis(ind) 116 130 ps=f_ps(ind) … … 125 139 out_u=f_out_u(ind) 126 140 !$OMP PARALLEL DEFAULT(SHARED) 127 ! CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du) 141 CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du) 142 !$OMP END PARALLEL 143 ENDDO 144 145 CASE(2) ! enstrophy-conserving 146 DO ind=1,ndomain 147 CALL swap_dimensions(ind) 148 CALL swap_geometry(ind) 149 phis=f_phis(ind) 150 ps=f_ps(ind) 151 dps=f_dps(ind) 152 theta_rhodz=f_theta_rhodz(ind) 153 dtheta_rhodz=f_dtheta_rhodz(ind) 154 rhodz=f_rhodz(ind) 155 p=f_p(ind) 156 qu=f_qu(ind) 157 u=f_u(ind) 158 du=f_du(ind) 159 out_u=f_out_u(ind) 160 !$OMP PARALLEL DEFAULT(SHARED) 128 161 CALL compute_pvort(ps, u, p,rhodz,qu) 129 162 CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du) … … 131 164 ENDDO 132 165 133 IF (mod(it,itau_out)==0 ) THEN134 PRINT *,'CALL write_output_fields'135 CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &136 f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p)137 END IF138 139 166 CASE DEFAULT 140 167 STOP 141 168 END SELECT 142 169 170 IF (mod(it,itau_out)==0 ) THEN 171 PRINT *,'CALL write_output_fields' 172 CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, & 173 f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p) 174 END IF 175 143 176 ! CALL check_mass_conservation(f_ps,f_dps) 144 177 145 178 END SUBROUTINE caldyn 146 147 179 148 180 SUBROUTINE compute_pvort(ps, u, p,rhodz,qu) 149 181 USE icosa … … 151 183 USE exner_mod 152 184 IMPLICIT NONE 153 REAL(rstd),INTENT(IN) :: u(iim*3*jjm,llm) 154 REAL(rstd),INTENT(IN) :: ps(iim*jjm) 155 REAL(rstd),INTENT(OUT) :: p(iim*jjm,llm+1) 156 REAL(rstd),INTENT(OUT) :: rhodz(iim*jjm,llm) 157 REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm) 158 159 INTEGER :: i,j,ij,l 160 ! REAL(rstd) :: ww,uu 161 ! REAL(rstd) :: delta 162 REAL(rstd) :: etav,hv 163 REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:) ! potential velocity 164 165 LOGICAL,SAVE :: first=.TRUE. 166 !$OMP THREADPRIVATE(first) 167 168 !$OMP BARRIER 169 !$OMP MASTER 170 ! IF (first) THEN 171 ALLOCATE(qv(2*iim*jjm,llm)) ! potential velocity 172 185 REAL(rstd),INTENT(IN) :: u(iim*3*jjm,llm) 186 REAL(rstd),INTENT(IN) :: ps(iim*jjm) 187 REAL(rstd),INTENT(OUT) :: p(iim*jjm,llm+1) 188 REAL(rstd),INTENT(OUT) :: rhodz(iim*jjm,llm) 189 REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm) 190 191 INTEGER :: i,j,ij,l 192 REAL(rstd) :: etav,hv 193 REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:) ! potential velocity 194 195 LOGICAL,SAVE :: first=.TRUE. 196 !$OMP THREADPRIVATE(first) 197 198 !$OMP BARRIER 199 !$OMP MASTER 200 ! IF (first) THEN 201 ALLOCATE(qv(2*iim*jjm,llm)) ! potential velocity 202 173 203 !!! Compute pressure 174 175 176 177 178 179 180 181 182 183 204 DO l = 1, llm+1 205 !$OMP DO 206 DO j=jj_begin-1,jj_end+1 207 DO i=ii_begin-1,ii_end+1 208 ij=(j-1)*iim+i 209 p(ij,l) = ap(l) + bp(l) * ps(ij) 210 ENDDO 211 ENDDO 212 ENDDO 213 184 214 !!! Compute mass 185 186 !$OMP DO187 188 189 190 191 192 193 194 215 DO l = 1, llm 216 !$OMP DO 217 DO j=jj_begin-1,jj_end+1 218 DO i=ii_begin-1,ii_end+1 219 ij=(j-1)*iim+i 220 rhodz(ij,l) = ( p(ij,l) - p(ij,l+1) )/g 221 ENDDO 222 ENDDO 223 ENDDO 224 195 225 !!! Compute shallow-water potential vorticity 196 226 DO l = 1,llm 197 227 !$OMP DO 198 228 DO j=jj_begin-1,jj_end+1 199 DO i=ii_begin-1,ii_end+1200 ij=(j-1)*iim+i201 229 DO i=ii_begin-1,ii_end+1 230 ij=(j-1)*iim+i 231 202 232 etav= 1./Av(ij+z_up)*( ne(ij,rup) * u(ij+u_rup,l) * de(ij+u_rup) & 203 233 + ne(ij+t_rup,left) * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left) & … … 236 266 !!$OMP BARRIER 237 267 !!$OMP MASTER 238 DEALLOCATE(theta) ! potential temperature239 DEALLOCATE(pk) ! Exner function240 DEALLOCATE(pks)241 DEALLOCATE(alpha)242 DEALLOCATE(beta)243 DEALLOCATE(phi) ! geopotential244 DEALLOCATE(Fe) ! mass flux245 DEALLOCATE(Ftheta) ! theta flux246 DEALLOCATE(convm) ! mass flux convergence247 DEALLOCATE(w) ! vertical velocity248 268 DEALLOCATE(qv) ! potential velocity 249 DEALLOCATE(berni) ! bernouilli term250 269 !!$OMP END MASTER 251 270 !!$OMP BARRIER 252 271 END SUBROUTINE compute_pvort 253 272 254 255 !----------- needed next : rhodz, p --------------!256 257 273 SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du) 258 274 USE icosa … … 284 300 REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:) ! mass flux convergence 285 301 REAL(rstd),ALLOCATABLE,SAVE :: w(:,:) ! vertical velocity 286 REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:) ! potential velocity287 302 REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:) ! Bernouilli function 288 303 … … 303 318 ALLOCATE(convm(iim*jjm,llm)) ! mass flux convergence 304 319 ALLOCATE(w(iim*jjm,llm)) ! vertical velocity 305 ALLOCATE(qv(2*iim*jjm,llm)) ! potential velocity306 320 ALLOCATE(berni(iim*jjm,llm)) ! bernouilli term 307 321 … … 394 408 395 409 !!! Compute potential vorticity (Coriolis) contribution to du 396 DO l=1,llm 397 !$OMP DO 398 DO j=jj_begin,jj_end 399 DO i=ii_begin,ii_end 400 ij=(j-1)*iim+i 401 402 du(ij+u_right,l) = qu(ij+u_right,l)/de(ij+u_right) * & 403 ( wee(ij+u_right,1,1)*Fe(ij+u_rup,l)+ & 404 wee(ij+u_right,2,1)*Fe(ij+u_lup,l)+ & 405 wee(ij+u_right,3,1)*Fe(ij+u_left,l)+ & 406 wee(ij+u_right,4,1)*Fe(ij+u_ldown,l)+ & 407 wee(ij+u_right,5,1)*Fe(ij+u_rdown,l)+ & 408 wee(ij+u_right,1,2)*Fe(ij+t_right+u_ldown,l)+ & 409 wee(ij+u_right,2,2)*Fe(ij+t_right+u_rdown,l)+ & 410 wee(ij+u_right,3,2)*Fe(ij+t_right+u_right,l)+ & 411 wee(ij+u_right,4,2)*Fe(ij+t_right+u_rup,l)+ & 412 wee(ij+u_right,5,2)*Fe(ij+t_right+u_lup,l) ) 413 414 415 du(ij+u_lup,l) = qu(ij+u_lup,l)/de(ij+u_lup) * & 416 ( wee(ij+u_lup,1,1)*Fe(ij+u_left,l)+ & 417 wee(ij+u_lup,2,1)*Fe(ij+u_ldown,l)+ & 418 wee(ij+u_lup,3,1)*Fe(ij+u_rdown,l)+ & 419 wee(ij+u_lup,4,1)*Fe(ij+u_right,l)+ & 420 wee(ij+u_lup,5,1)*Fe(ij+u_rup,l)+ & 421 wee(ij+u_lup,1,2)*Fe(ij+t_lup+u_right,l)+ & 422 wee(ij+u_lup,2,2)*Fe(ij+t_lup+u_rup,l)+ & 423 wee(ij+u_lup,3,2)*Fe(ij+t_lup+u_lup,l)+ & 424 wee(ij+u_lup,4,2)*Fe(ij+t_lup+u_left,l)+ & 425 wee(ij+u_lup,5,2)*Fe(ij+t_lup+u_ldown,l) ) 426 427 428 du(ij+u_ldown,l) = qu(ij+u_ldown,l)/de(ij+u_ldown) * & 429 ( wee(ij+u_ldown,1,1)*Fe(ij+u_rdown,l)+ & 430 wee(ij+u_ldown,2,1)*Fe(ij+u_right,l)+ & 431 wee(ij+u_ldown,3,1)*Fe(ij+u_rup,l)+ & 432 wee(ij+u_ldown,4,1)*Fe(ij+u_lup,l)+ & 433 wee(ij+u_ldown,5,1)*Fe(ij+u_left,l)+ & 434 wee(ij+u_ldown,1,2)*Fe(ij+t_ldown+u_lup,l)+ & 435 wee(ij+u_ldown,2,2)*Fe(ij+t_ldown+u_left,l)+ & 436 wee(ij+u_ldown,3,2)*Fe(ij+t_ldown+u_ldown,l)+ & 437 wee(ij+u_ldown,4,2)*Fe(ij+t_ldown+u_rdown,l)+ & 438 wee(ij+u_ldown,5,2)*Fe(ij+t_ldown+u_right,l) ) 439 410 411 SELECT CASE(caldyn_conserv) 412 CASE(1) ! energy-conserving TRiSK 413 414 DO l=1,llm 415 !$OMP DO 416 DO j=jj_begin,jj_end 417 DO i=ii_begin,ii_end 418 ij=(j-1)*iim+i 419 420 uu = wee(ij+u_right,1,1)*Fe(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+ & 421 wee(ij+u_right,2,1)*Fe(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+ & 422 wee(ij+u_right,3,1)*Fe(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+ & 423 wee(ij+u_right,4,1)*Fe(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+ & 424 wee(ij+u_right,5,1)*Fe(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+ & 425 wee(ij+u_right,1,2)*Fe(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+ & 426 wee(ij+u_right,2,2)*Fe(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+ & 427 wee(ij+u_right,3,2)*Fe(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+ & 428 wee(ij+u_right,4,2)*Fe(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+ & 429 wee(ij+u_right,5,2)*Fe(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l)) 430 du(ij+u_right,l) = .5*uu/de(ij+u_right) 431 432 uu = wee(ij+u_lup,1,1)*Fe(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) + & 433 wee(ij+u_lup,2,1)*Fe(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) + & 434 wee(ij+u_lup,3,1)*Fe(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) + & 435 wee(ij+u_lup,4,1)*Fe(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) + & 436 wee(ij+u_lup,5,1)*Fe(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) + & 437 wee(ij+u_lup,1,2)*Fe(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) + & 438 wee(ij+u_lup,2,2)*Fe(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) + & 439 wee(ij+u_lup,3,2)*Fe(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) + & 440 wee(ij+u_lup,4,2)*Fe(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) + & 441 wee(ij+u_lup,5,2)*Fe(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l)) 442 du(ij+u_lup,l) = .5*uu/de(ij+u_lup) 443 444 445 uu = wee(ij+u_ldown,1,1)*Fe(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) + & 446 wee(ij+u_ldown,2,1)*Fe(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) + & 447 wee(ij+u_ldown,3,1)*Fe(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) + & 448 wee(ij+u_ldown,4,1)*Fe(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) + & 449 wee(ij+u_ldown,5,1)*Fe(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) + & 450 wee(ij+u_ldown,1,2)*Fe(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) + & 451 wee(ij+u_ldown,2,2)*Fe(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) + & 452 wee(ij+u_ldown,3,2)*Fe(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) + & 453 wee(ij+u_ldown,4,2)*Fe(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) + & 454 wee(ij+u_ldown,5,2)*Fe(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l)) 455 du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown) 456 457 ENDDO 458 ENDDO 459 ENDDO 460 461 CASE(2) ! enstrophy-conserving TRiSK 462 463 DO l=1,llm 464 !$OMP DO 465 DO j=jj_begin,jj_end 466 DO i=ii_begin,ii_end 467 ij=(j-1)*iim+i 468 469 uu = wee(ij+u_right,1,1)*Fe(ij+u_rup,l)+ & 470 wee(ij+u_right,2,1)*Fe(ij+u_lup,l)+ & 471 wee(ij+u_right,3,1)*Fe(ij+u_left,l)+ & 472 wee(ij+u_right,4,1)*Fe(ij+u_ldown,l)+ & 473 wee(ij+u_right,5,1)*Fe(ij+u_rdown,l)+ & 474 wee(ij+u_right,1,2)*Fe(ij+t_right+u_ldown,l)+ & 475 wee(ij+u_right,2,2)*Fe(ij+t_right+u_rdown,l)+ & 476 wee(ij+u_right,3,2)*Fe(ij+t_right+u_right,l)+ & 477 wee(ij+u_right,4,2)*Fe(ij+t_right+u_rup,l)+ & 478 wee(ij+u_right,5,2)*Fe(ij+t_right+u_lup,l) 479 du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right) 480 481 482 uu = wee(ij+u_lup,1,1)*Fe(ij+u_left,l)+ & 483 wee(ij+u_lup,2,1)*Fe(ij+u_ldown,l)+ & 484 wee(ij+u_lup,3,1)*Fe(ij+u_rdown,l)+ & 485 wee(ij+u_lup,4,1)*Fe(ij+u_right,l)+ & 486 wee(ij+u_lup,5,1)*Fe(ij+u_rup,l)+ & 487 wee(ij+u_lup,1,2)*Fe(ij+t_lup+u_right,l)+ & 488 wee(ij+u_lup,2,2)*Fe(ij+t_lup+u_rup,l)+ & 489 wee(ij+u_lup,3,2)*Fe(ij+t_lup+u_lup,l)+ & 490 wee(ij+u_lup,4,2)*Fe(ij+t_lup+u_left,l)+ & 491 wee(ij+u_lup,5,2)*Fe(ij+t_lup+u_ldown,l) 492 du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup) 493 494 uu = wee(ij+u_ldown,1,1)*Fe(ij+u_rdown,l)+ & 495 wee(ij+u_ldown,2,1)*Fe(ij+u_right,l)+ & 496 wee(ij+u_ldown,3,1)*Fe(ij+u_rup,l)+ & 497 wee(ij+u_ldown,4,1)*Fe(ij+u_lup,l)+ & 498 wee(ij+u_ldown,5,1)*Fe(ij+u_left,l)+ & 499 wee(ij+u_ldown,1,2)*Fe(ij+t_ldown+u_lup,l)+ & 500 wee(ij+u_ldown,2,2)*Fe(ij+t_ldown+u_left,l)+ & 501 wee(ij+u_ldown,3,2)*Fe(ij+t_ldown+u_ldown,l)+ & 502 wee(ij+u_ldown,4,2)*Fe(ij+t_ldown+u_rdown,l)+ & 503 wee(ij+u_ldown,5,2)*Fe(ij+t_ldown+u_right,l) 504 du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown) 440 505 441 ENDDO 442 ENDDO 443 ENDDO 506 ENDDO 507 ENDDO 508 ENDDO 509 510 CASE DEFAULT 511 STOP 512 END SELECT 444 513 445 514 !!! Compute Exner function … … 567 636 DEALLOCATE(convm) ! mass flux convergence 568 637 DEALLOCATE(w) ! vertical velocity 569 DEALLOCATE(qv) ! potential velocity570 638 DEALLOCATE(berni) ! bernouilli term 571 639 !!$OMP END MASTER
Note: See TracChangeset
for help on using the changeset viewer.
