![(please configure the [header_logo] section in trac.ini)](/dynamico/chrome/site/your_project_logo.png){kind=logo}
Changeset 162 for codes/icosagcm/trunk/src/caldyn_gcm.f90
- Timestamp:
- 06/27/13 18:37:27 (11 years ago)
- File:
-
- 1 edited
-
codes/icosagcm/trunk/src/caldyn_gcm.f90 (modified) (26 diffs)
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/caldyn_gcm.f90
r159 r162 5 5 6 6 INTEGER, PARAMETER :: energy=1, enstrophy=2 7 TYPE(t_field),POINTER :: f_out_u(:), f_qu(:) 7 TYPE(t_field),POINTER :: f_out_u(:), f_qu(:), f_qv(:) 8 8 REAL(rstd),POINTER :: out_u(:,:), p(:,:), qu(:,:) 9 9 … … 15 15 TYPE(t_field),POINTER :: f_pk(:) 16 16 TYPE(t_field),POINTER :: f_geopot(:) 17 TYPE(t_field),POINTER :: f_divm(:)18 17 TYPE(t_field),POINTER :: f_wwuu(:) 19 18 20 19 INTEGER :: caldyn_conserv 21 20 22 TYPE(t_message) :: req_ps, req_ theta_rhodz, req_u, req_qu23 24 PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields, req_ps,&25 caldyn_eta, eta_mass, eta_lag21 TYPE(t_message) :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu 22 23 PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields, & 24 req_ps, req_mass, caldyn_eta, eta_mass, eta_lag 26 25 27 26 CONTAINS … … 58 57 CALL allocate_field(f_out_u,field_u,type_real,llm) 59 58 CALL allocate_field(f_qu,field_u,type_real,llm) 59 CALL allocate_field(f_qv,field_z,type_real,llm) 60 60 61 61 CALL allocate_field(f_buf_i, field_t,type_real,llm) … … 70 70 CALL allocate_field(f_pk, field_t,type_real,llm, name='pk') 71 71 CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot') ! geopotential 72 CALL allocate_field(f_divm, field_t,type_real,llm, name='divm') ! mass flux divergence73 72 CALL allocate_field(f_wwuu, field_u,type_real,llm+1,name='wwuu') 74 73 … … 121 120 122 121 SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, & 123 f_hflux, f_wflux, f_dps, f_d theta_rhodz, f_du)122 f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du) 124 123 USE icosa 124 USE disvert_mod 125 125 USE vorticity_mod 126 126 USE kinetic_mod … … 139 139 TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:) 140 140 TYPE(t_field),POINTER :: f_dps(:) 141 TYPE(t_field),POINTER :: f_dmass(:) 141 142 TYPE(t_field),POINTER :: f_dtheta_rhodz(:) 142 143 TYPE(t_field),POINTER :: f_du(:) … … 146 147 REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:) 147 148 REAL(rstd),POINTER :: qu(:,:) 149 REAL(rstd),POINTER :: qv(:,:) 148 150 149 151 ! temporary shared variable … … 151 153 REAL(rstd),POINTER :: pk(:,:) 152 154 REAL(rstd),POINTER :: geopot(:,:) 153 REAL(rstd),POINTER :: divm(:,:)155 REAL(rstd),POINTER :: convm(:,:) 154 156 REAL(rstd),POINTER :: wwuu(:,:) 155 157 … … 158 160 !$OMP THREADPRIVATE(first) 159 161 160 IF (first) THEN 162 ! MPI messages need to be sent at first call to caldyn 163 ! This is needed only once : the next ones will be sent by timeloop 164 IF (first) THEN 161 165 first=.FALSE. 162 CALL init_message(f_ps,req_i1,req_ps) 166 IF(caldyn_eta==eta_mass) THEN 167 CALL init_message(f_ps,req_i1,req_ps) 168 ELSE 169 CALL init_message(f_mass,req_i1,req_mass) 170 END IF 163 171 CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz) 164 172 CALL init_message(f_u,req_e1_vect,req_u) 165 173 CALL init_message(f_qu,req_e1_scal,req_qu) 166 CALL send_message(f_ps,req_ps) 174 IF(caldyn_eta==eta_mass) THEN 175 CALL send_message(f_ps,req_ps) 176 ELSE 177 CALL send_message(f_mass,req_mass) 178 END IF 167 179 ENDIF 168 180 … … 183 195 theta = f_theta(ind) 184 196 qu=f_qu(ind) 185 CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu) 197 qv=f_qv(ind) 198 CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv) 186 199 ENDDO 187 200 … … 201 214 CALL compute_geopot(ps,mass,theta, pk,geopot) 202 215 hflux=f_hflux(ind) 203 divm = f_divm(ind)216 convm = f_dmass(ind) 204 217 dtheta_rhodz=f_dtheta_rhodz(ind) 205 218 du=f_du(ind) 206 CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du) 207 wflux=f_wflux(ind) 208 wwuu=f_wwuu(ind) 209 dps=f_dps(ind) 210 CALL compute_caldyn_vert(u,theta,mass,divm, wflux,wwuu, dps, dtheta_rhodz, du) 219 CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du) 220 IF(caldyn_eta==eta_mass) THEN 221 wflux=f_wflux(ind) 222 wwuu=f_wwuu(ind) 223 dps=f_dps(ind) 224 CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du) 225 END IF 211 226 ENDDO 212 227 … … 221 236 theta = f_theta(ind) 222 237 qu=f_qu(ind) 223 CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu )238 CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv) 224 239 pk = f_pk(ind) 225 240 geopot = f_geopot(ind) 226 241 CALL compute_geopot(ps,mass,theta, pk,geopot) 227 242 hflux=f_hflux(ind) 228 divm = f_divm(ind)243 convm = f_dmass(ind) 229 244 dtheta_rhodz=f_dtheta_rhodz(ind) 230 245 du=f_du(ind) 231 CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du) 232 wflux=f_wflux(ind) 233 wwuu=f_wwuu(ind) 234 dps=f_dps(ind) 235 CALL compute_caldyn_vert(u,theta,mass,divm, wflux,wwuu, dps, dtheta_rhodz, du) 246 CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du) 247 IF(caldyn_eta==eta_mass) THEN 248 wflux=f_wflux(ind) 249 wwuu=f_wwuu(ind) 250 dps=f_dps(ind) 251 CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du) 252 END IF 236 253 ENDDO 237 254 … … 252 269 CALL writefield("dps",f_dps) 253 270 CALL writefield("mass",f_mass) 271 CALL writefield("dmass",f_dmass) 254 272 CALL writefield("vort",f_qu) 255 273 CALL writefield("theta",f_theta) 274 CALL writefield("exner",f_pk) 275 CALL writefield("pv",f_qv) 256 276 257 277 END IF … … 263 283 END SUBROUTINE caldyn 264 284 265 SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu )285 SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv) 266 286 USE icosa 267 287 USE disvert_mod … … 273 293 REAL(rstd),INTENT(IN) :: ps(iim*jjm) 274 294 REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm) 275 REAL(rstd),INTENT( OUT) :: rhodz(iim*jjm,llm)295 REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm) 276 296 REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm) 277 297 REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm) 298 REAL(rstd),INTENT(OUT) :: qv(iim*2*jjm,llm) 278 299 279 300 INTEGER :: i,j,ij,l 280 REAL(rstd) :: etav,hv 281 REAL(rstd) :: qv(2*iim*jjm,llm) ! potential velocity301 REAL(rstd) :: etav,hv, m 302 ! REAL(rstd) :: qv(2*iim*jjm,llm) ! potential velocity 282 303 283 304 284 305 CALL trace_start("compute_pvort") 285 306 286 CALL wait_message(req_ps) 307 IF(caldyn_eta==eta_mass) THEN 308 CALL wait_message(req_ps) 309 ELSE 310 CALL wait_message(req_mass) 311 END IF 287 312 CALL wait_message(req_theta_rhodz) 288 313 289 !!! Compute mass & theta 290 DO l = ll_begin,ll_end 291 CALL test_message(req_u) 292 DO j=jj_begin-1,jj_end+1 293 DO i=ii_begin-1,ii_end+1 294 ij=(j-1)*iim+i 295 rhodz(ij,l) = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g 296 theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 314 IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta 315 DO l = ll_begin,ll_end 316 CALL test_message(req_u) 317 DO j=jj_begin-1,jj_end+1 318 DO i=ii_begin-1,ii_end+1 319 ij=(j-1)*iim+i 320 m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 321 ! m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g 322 rhodz(ij,l) = m 323 theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 324 ENDDO 297 325 ENDDO 298 326 ENDDO 299 ENDDO 327 ELSE ! Compute only theta 328 DO l = ll_begin,ll_end 329 CALL test_message(req_u) 330 DO j=jj_begin-1,jj_end+1 331 DO i=ii_begin-1,ii_end+1 332 ij=(j-1)*iim+i 333 theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 334 ENDDO 335 ENDDO 336 ENDDO 337 END IF 300 338 301 339 CALL wait_message(req_u) … … 377 415 pk(ij,l) = exner_ik 378 416 ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 379 geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 417 geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 380 418 ENDDO 381 419 ENDDO … … 413 451 END DO 414 452 END DO 415 416 453 DO l = 1,llm 417 454 !$OMP DO SCHEDULE(STATIC) … … 434 471 END SUBROUTINE compute_geopot 435 472 436 SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux, divm, dtheta_rhodz, du)473 SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,convm, dtheta_rhodz, du) 437 474 USE icosa 438 475 USE disvert_mod … … 449 486 450 487 REAL(rstd),INTENT(OUT) :: hflux(iim*3*jjm,llm) ! hflux in kg/s 451 REAL(rstd),INTENT(OUT) :: divm(iim*jjm,llm) ! mass flux divergence488 REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm) ! mass flux convergence 452 489 REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm) 453 490 REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm) … … 483 520 DO i=ii_begin,ii_end 484 521 ij=(j-1)*iim+i 485 ! divm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21486 divm(ij,l)=1./Ai(ij)*(ne_right*hflux(ij+u_right,l) + &522 ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 523 convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l) + & 487 524 ne_rup*hflux(ij+u_rup,l) + & 488 525 ne_lup*hflux(ij+u_lup,l) + & … … 659 696 END SUBROUTINE compute_caldyn_horiz 660 697 661 SUBROUTINE compute_caldyn_vert(u,theta,rhodz, divm, wflux,wwuu, dps,dtheta_rhodz,du)698 SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du) 662 699 USE icosa 663 700 USE disvert_mod … … 669 706 REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm) 670 707 REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 671 REAL(rstd),INTENT(INOUT) :: divm(iim*jjm,llm) ! mass flux divergence708 REAL(rstd),INTENT(INOUT) :: convm(iim*jjm,llm) ! mass flux convergence 672 709 673 710 REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s) … … 684 721 685 722 !$OMP BARRIER 686 !!! cumulate mass flux divergence from top to bottom723 !!! cumulate mass flux convergence from top to bottom 687 724 DO l = llm-1, 1, -1 688 725 IF (caldyn_conserv==energy) CALL test_message(req_qu) … … 691 728 DO i=ii_begin,ii_end 692 729 ij=(j-1)*iim+i 693 divm(ij,l) = divm(ij,l) + divm(ij,l+1)730 convm(ij,l) = convm(ij,l) + convm(ij,l+1) 694 731 ENDDO 695 732 ENDDO 696 733 ENDDO 697 734 698 ! IMPLICIT FLUSH on divm735 ! IMPLICIT FLUSH on convm 699 736 !!!!!!!!!!!!!!!!!!!!!!!!! 700 737 … … 705 742 ij=(j-1)*iim+i 706 743 ! dps/dt = -int(div flux)dz 707 dps(ij) =-divm(ij,1) * g744 dps(ij) = convm(ij,1) * g 708 745 ENDDO 709 746 ENDDO … … 718 755 ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt 719 756 ! => w>0 for upward transport 720 wflux( ij, l ) = divm( ij, l ) - bp(l) * divm( ij, 1 )757 wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l ) 721 758 ENDDO 722 759 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
