Changeset 159
- Timestamp:
- 06/25/13 09:32:16 (11 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/caldyn.f90
r157 r159 31 31 END SUBROUTINE init_caldyn 32 32 33 SUBROUTINE caldyn(write_out,f_phis, f_ps, f_ theta_rhodz, f_u, f_q, &33 SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, & 34 34 f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) 35 35 USE icosa … … 40 40 TYPE(t_field),POINTER :: f_phis(:) 41 41 TYPE(t_field),POINTER :: f_ps(:) 42 TYPE(t_field),POINTER :: f_mass(:) 42 43 TYPE(t_field),POINTER :: f_theta_rhodz(:) 43 44 TYPE(t_field),POINTER :: f_u(:) … … 51 52 SELECT CASE (TRIM(caldyn_type)) 52 53 CASE('gcm') 53 CALL caldyn_gcm(write_out,f_phis, f_ps, f_ theta_rhodz, f_u, f_q, &54 CALL caldyn_gcm(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, & 54 55 f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) 55 56 CASE('adv') -
codes/icosagcm/trunk/src/caldyn_gcm.f90
r157 r159 5 5 6 6 INTEGER, PARAMETER :: energy=1, enstrophy=2 7 TYPE(t_field),POINTER :: f_out_u(:), f_ rhodz(:), f_qu(:)8 REAL(rstd),POINTER :: out_u(:,:), p(:,:), rhodz(:,:),qu(:,:)7 TYPE(t_field),POINTER :: f_out_u(:), f_qu(:) 8 REAL(rstd),POINTER :: out_u(:,:), p(:,:), qu(:,:) 9 9 10 10 TYPE(t_field),POINTER :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:) … … 18 18 TYPE(t_field),POINTER :: f_wwuu(:) 19 19 20 INTEGER :: caldyn_ hydrostat, caldyn_conserv20 INTEGER :: caldyn_conserv 21 21 22 22 TYPE(t_message) :: req_ps, req_theta_rhodz, req_u, req_qu 23 23 24 PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields,req_ps 24 PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields,req_ps, & 25 caldyn_eta, eta_mass, eta_lag 25 26 26 27 CONTAINS … … 41 42 caldyn_conserv=enstrophy 42 43 CASE DEFAULT 43 IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', TRIM(def),'> options are <energy>, <enstrophy>' 44 IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', & 45 TRIM(def),'> options are <energy>, <enstrophy>' 44 46 STOP 45 47 END SELECT … … 55 57 56 58 CALL allocate_field(f_out_u,field_u,type_real,llm) 57 CALL allocate_field(f_rhodz,field_t,type_real,llm)58 59 CALL allocate_field(f_qu,field_u,type_real,llm) 59 60 60 CALL allocate_field(f_buf_i, field_t,type_real,llm)61 CALL allocate_field(f_buf_p, field_t,type_real,llm+1)62 CALL allocate_field(f_buf_u3d, field_t,type_real,3,llm) ! 3D vel at cell centers61 CALL allocate_field(f_buf_i, field_t,type_real,llm) 62 CALL allocate_field(f_buf_p, field_t,type_real,llm+1) 63 CALL allocate_field(f_buf_u3d, field_t,type_real,3,llm) ! 3D vel at cell centers 63 64 CALL allocate_field(f_buf_ulon,field_t,type_real,llm) 64 65 CALL allocate_field(f_buf_ulat,field_t,type_real,llm) 65 CALL allocate_field(f_buf_v, field_z,type_real,llm)66 CALL allocate_field(f_buf_s, field_t,type_real)67 68 CALL allocate_field(f_theta, field_t,type_real,llm)! potential temperature69 CALL allocate_field(f_pk, field_t,type_real,llm)70 CALL allocate_field(f_geopot,field_t,type_real,llm+1 ) ! geopotential71 CALL allocate_field(f_divm, field_t,type_real,llm)! mass flux divergence72 CALL allocate_field(f_wwuu, field_u,type_real,llm+1)66 CALL allocate_field(f_buf_v, field_z,type_real,llm) 67 CALL allocate_field(f_buf_s, field_t,type_real) 68 69 CALL allocate_field(f_theta, field_t,type_real,llm, name='theta') ! potential temperature 70 CALL allocate_field(f_pk, field_t,type_real,llm, name='pk') 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 divergence 73 CALL allocate_field(f_wwuu, field_u,type_real,llm+1,name='wwuu') 73 74 74 75 END SUBROUTINE allocate_caldyn … … 119 120 END SUBROUTINE caldyn_BC 120 121 121 SUBROUTINE caldyn(write_out,f_phis, f_ps, f_ theta_rhodz, f_u, f_q, &122 SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, & 122 123 f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) 123 124 USE icosa … … 132 133 TYPE(t_field),POINTER :: f_phis(:) 133 134 TYPE(t_field),POINTER :: f_ps(:) 135 TYPE(t_field),POINTER :: f_mass(:) 134 136 TYPE(t_field),POINTER :: f_theta_rhodz(:) 135 137 TYPE(t_field),POINTER :: f_u(:) … … 140 142 TYPE(t_field),POINTER :: f_du(:) 141 143 142 REAL(rstd),POINTER :: p his(:), ps(:), dps(:)143 REAL(rstd),POINTER :: theta_rhodz(:,:), dtheta_rhodz(:,:)144 REAL(rstd),POINTER :: ps(:), dps(:) 145 REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:), dtheta_rhodz(:,:) 144 146 REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:) 145 REAL(rstd),POINTER :: rhodz(:,:),qu(:,:)147 REAL(rstd),POINTER :: qu(:,:) 146 148 147 149 ! temporary shared variable … … 178 180 u=f_u(ind) 179 181 theta_rhodz = f_theta_rhodz(ind) 180 rhodz=f_rhodz(ind)182 mass=f_mass(ind) 181 183 theta = f_theta(ind) 182 184 qu=f_qu(ind) 183 CALL compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu)185 CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu) 184 186 ENDDO 185 187 … … 192 194 u=f_u(ind) 193 195 theta_rhodz=f_theta_rhodz(ind) 194 rhodz=f_rhodz(ind)196 mass=f_mass(ind) 195 197 theta = f_theta(ind) 196 198 qu=f_qu(ind) 197 phis=f_phis(ind)198 199 pk = f_pk(ind) 199 200 geopot = f_geopot(ind) 200 CALL compute_geopot(ps, phis,rhodz,theta, pk,geopot)201 CALL compute_geopot(ps,mass,theta, pk,geopot) 201 202 hflux=f_hflux(ind) 202 203 divm = f_divm(ind) 203 204 dtheta_rhodz=f_dtheta_rhodz(ind) 204 205 du=f_du(ind) 205 CALL compute_caldyn_horiz(u, rhodz,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du)206 CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du) 206 207 wflux=f_wflux(ind) 207 208 wwuu=f_wwuu(ind) 208 209 dps=f_dps(ind) 209 CALL compute_caldyn_vert(u,theta, rhodz,divm, wflux,wwuu, dps, dtheta_rhodz, du)210 CALL compute_caldyn_vert(u,theta,mass,divm, wflux,wwuu, dps, dtheta_rhodz, du) 210 211 ENDDO 211 212 … … 217 218 u=f_u(ind) 218 219 theta_rhodz=f_theta_rhodz(ind) 219 rhodz=f_rhodz(ind)220 mass=f_mass(ind) 220 221 theta = f_theta(ind) 221 222 qu=f_qu(ind) 222 CALL compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu) 223 phis=f_phis(ind) 223 CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu) 224 224 pk = f_pk(ind) 225 225 geopot = f_geopot(ind) 226 CALL compute_geopot(ps, phis,rhodz,theta, pk,geopot)226 CALL compute_geopot(ps,mass,theta, pk,geopot) 227 227 hflux=f_hflux(ind) 228 228 divm = f_divm(ind) 229 229 dtheta_rhodz=f_dtheta_rhodz(ind) 230 230 du=f_du(ind) 231 CALL compute_caldyn_horiz(u, rhodz,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du)231 CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du) 232 232 wflux=f_wflux(ind) 233 233 wwuu=f_wwuu(ind) 234 234 dps=f_dps(ind) 235 CALL compute_caldyn_vert(u,theta, rhodz,divm, wflux,wwuu, dps, dtheta_rhodz, du)235 CALL compute_caldyn_vert(u,theta,mass,divm, wflux,wwuu, dps, dtheta_rhodz, du) 236 236 ENDDO 237 237 … … 248 248 ! CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, & 249 249 ! f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p) 250 250 251 CALL writefield("ps",f_ps) 251 252 CALL writefield("dps",f_dps) 253 CALL writefield("mass",f_mass) 252 254 CALL writefield("vort",f_qu) 253 255 CALL writefield("theta",f_theta) … … 344 346 END SUBROUTINE compute_pvort 345 347 346 SUBROUTINE compute_geopot(ps, phis,rhodz,theta,pk,geopot)348 SUBROUTINE compute_geopot(ps,rhodz,theta,pk,geopot) 347 349 USE icosa 348 350 USE disvert_mod … … 351 353 USE omp_para 352 354 IMPLICIT NONE 353 REAL(rstd),INTENT(IN) :: phis(iim*jjm) 354 REAL(rstd),INTENT(IN) :: ps(iim*jjm) 355 REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 356 REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm) ! potential temperature 357 REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function 358 REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential 355 REAL(rstd),INTENT(INOUT) :: ps(iim*jjm) 356 REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 357 REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm) ! potential temperature 358 REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function 359 REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential 359 360 360 361 INTEGER :: i,j,ij,l … … 363 364 CALL trace_start("compute_geopot") 364 365 366 IF(caldyn_eta==eta_mass) THEN 367 365 368 !!! Compute exner function and geopotential 366 DO l = 1,llm 367 !$OMP DO SCHEDULE(STATIC) 369 DO l = 1,llm 370 !$OMP DO SCHEDULE(STATIC) 371 DO j=jj_begin-1,jj_end+1 372 DO i=ii_begin-1,ii_end+1 373 ij=(j-1)*iim+i 374 p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later 375 ! p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j)) 376 exner_ik = cpp * (p_ik/preff) ** kappa 377 pk(ij,l) = exner_ik 378 ! 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 380 ENDDO 381 ENDDO 382 ENDDO 383 384 ELSE 385 386 ! We are using a Lagrangian vertical coordinate 387 ! Pressure must be computed first top-down (temporarily stored in pk) 388 ! Then Exner pressure and geopotential are computed bottom-up 389 ! Notice that the computation below should work also when caldyn_eta=eta_mass 390 391 ! uppermost layer 368 392 DO j=jj_begin-1,jj_end+1 369 393 DO i=ii_begin-1,ii_end+1 370 394 ij=(j-1)*iim+i 371 p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later 372 ! p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j)) 373 exner_ik = cpp * (p_ik/preff) ** kappa 374 pk(ij,l) = exner_ik 375 ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 376 geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 395 pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm) 396 END DO 397 END DO 398 ! other layers 399 DO l = llm-1, 1, -1 400 !$OMP DO SCHEDULE(STATIC) 401 DO j=jj_begin-1,jj_end+1 402 DO i=ii_begin-1,ii_end+1 403 ij=(j-1)*iim+i 404 pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1)) 405 END DO 406 END DO 407 END DO 408 ! surface pressure (for diagnostics) 409 DO j=jj_begin-1,jj_end+1 410 DO i=ii_begin-1,ii_end+1 411 ij=(j-1)*iim+i 412 ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1) 413 END DO 414 END DO 415 416 DO l = 1,llm 417 !$OMP DO SCHEDULE(STATIC) 418 DO j=jj_begin-1,jj_end+1 419 DO i=ii_begin-1,ii_end+1 420 ij=(j-1)*iim+i 421 p_ik = pk(ij,l) 422 exner_ik = cpp * (p_ik/preff) ** kappa 423 ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 424 geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 425 pk(ij,l) = exner_ik 426 ENDDO 377 427 ENDDO 378 428 ENDDO 379 ENDDO 429 430 END IF 380 431 381 432 CALL trace_end("compute_geopot") … … 798 849 799 850 ! Temperature 800 CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; 801 851 ! CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; ! FIXME 852 802 853 CALL getin('physics',physics_type) 803 854 IF (TRIM(physics_type)=='dcmip') THEN … … 828 879 ENDIF 829 880 830 ! geopotential 881 ! geopotential ! FIXME 831 882 CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i) 832 883 CALL writefield("p",f_buf_p) 833 CALL writefield("phi",f_ buf_i)884 CALL writefield("phi",f_geopot) ! geopotential 834 885 CALL writefield("theta",f_buf1_i) ! potential temperature 835 CALL writefield("pk",f_buf2_i) ! Exner pressure 836 886 CALL writefield("pk",f_buf2_i) ! Exner pressure 837 887 838 888 END SUBROUTINE write_output_fields -
codes/icosagcm/trunk/src/disvert.f90
r156 r159 5 5 REAL(rstd), SAVE, POINTER :: presnivs(:) 6 6 REAL(rstd), SAVE, POINTER :: mass_al(:), mass_bl(:), mass_ak(:), mass_bk(:), mass_dak(:), mass_dbk(:) 7 7 8 REAL(rstd) :: ptop ! pressure at top of atmosphere l=llm+1 9 ! vertical coordinate : Lagrangian or mass-based 10 INTEGER, SAVE :: caldyn_eta 11 INTEGER, PARAMETER :: eta_mass=1, eta_lag=2 8 12 9 13 CONTAINS … … 91 95 END SUBROUTINE init_disvert 92 96 97 SUBROUTINE compute_rhodz(comp, ps, rhodz) 98 USE icosa 99 USE omp_para 100 LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check 101 REAL(rstd), INTENT(IN) :: ps(iim*jjm) 102 REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm) 103 REAL(rstd) :: m, err 104 INTEGER :: l,i,j,ij,dd 105 err=0. 106 107 dd=0 108 109 DO l = ll_begin, ll_end 110 DO j=jj_begin-dd,jj_end+dd 111 DO i=ii_begin-dd,ii_end+dd 112 ij=(j-1)*iim+i 113 m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 114 IF(comp) THEN 115 rhodz(ij,l) = m 116 ELSE 117 err = MAX(err,abs(m-rhodz(ij,l))) 118 END IF 119 ENDDO 120 ENDDO 121 ENDDO 122 123 IF(.NOT. comp) THEN 124 IF(err>1e-10) THEN 125 PRINT *, 'Discrepancy between ps and rhodz detected', err 126 STOP 127 ELSE 128 ! PRINT *, 'No discrepancy between ps and rhodz detected' 129 END IF 130 END IF 131 132 END SUBROUTINE compute_rhodz 133 134 93 135 SUBROUTINE write_apbp 94 136 USE icosa -
codes/icosagcm/trunk/src/etat0.f90
r149 r159 4 4 CONTAINS 5 5 6 SUBROUTINE etat0(f_ps,f_ phis,f_theta_rhodz,f_u, f_q)6 SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q) 7 7 USE icosa 8 USE disvert_mod 8 9 USE etat0_jablonowsky06_mod, ONLY : etat0_jablonowsky06=>etat0 9 10 USE etat0_academic_mod, ONLY : etat0_academic=>etat0 … … 19 20 IMPLICIT NONE 20 21 TYPE(t_field),POINTER :: f_ps(:) 22 TYPE(t_field),POINTER :: f_mass(:) 21 23 TYPE(t_field),POINTER :: f_phis(:) 22 24 TYPE(t_field),POINTER :: f_theta_rhodz(:) 23 25 TYPE(t_field),POINTER :: f_u(:) 24 26 TYPE(t_field),POINTER :: f_q(:) 25 27 REAL(rstd),POINTER :: ps(:), mass(:,:) 28 LOGICAL :: init_mass 29 INTEGER :: ind,i,j,ij,l 30 31 ! most etat0 routines set ps and not mass 32 ! in that case and if caldyn_eta == eta_lag 33 ! the initial distribution of mass is taken to be the same 34 ! as what the mass coordinate would dictate 35 ! however if etat0_XXX defines mass then the flag init_mass must be set to .FALSE. 36 ! otherwise mass will be overwritten 37 init_mass = (caldyn_eta == eta_lag) 38 26 39 etat0_type='jablonowsky06' 27 40 CALL getin("etat0",etat0_type) … … 56 69 STOP 57 70 END SELECT 58 71 72 IF(init_mass) THEN ! initialize mass distribution using ps 73 !$OMP BARRIER 74 DO ind=1,ndomain 75 CALL swap_dimensions(ind) 76 CALL swap_geometry(ind) 77 mass=f_mass(ind); ps=f_ps(ind) 78 CALL compute_rhodz(.TRUE., ps, mass) 79 END DO 80 END IF 81 59 82 END SUBROUTINE etat0 60 83 61 84 END MODULE etat0_mod 85 -
codes/icosagcm/trunk/src/field.f90
r138 r159 63 63 field(ind)%name = name 64 64 ELSE 65 field(ind)%name = '(un kown)'65 field(ind)%name = '(undefined)' 66 66 END IF 67 67 … … 232 232 233 233 IF (field%ndim/=2 .OR. field%data_type/=type_real) THEN 234 PRINT *, 'get_val_r2d : bad pointer assign ationwith ' // TRIM(field%name)234 PRINT *, 'get_val_r2d : bad pointer assignment with ' // TRIM(field%name) 235 235 STOP 236 236 END IF … … 244 244 245 245 IF (field%ndim/=3 .OR. field%data_type/=type_real) THEN 246 PRINT *, 'get_val_r3d : bad pointer assign ationwith ' // TRIM(field%name)246 PRINT *, 'get_val_r3d : bad pointer assignment with ' // TRIM(field%name) 247 247 STOP 248 ! CALL ABORT 248 249 END IF 249 250 field_pt=>field%rval3d … … 256 257 257 258 IF (field%ndim/=4 .OR. field%data_type/=type_real) THEN 258 PRINT *, 'get_val_r4d : bad pointer assign ationwith ' // TRIM(field%name)259 PRINT *, 'get_val_r4d : bad pointer assignment with ' // TRIM(field%name) 259 260 STOP 260 261 END IF … … 268 269 TYPE(t_field),INTENT(IN) :: field 269 270 270 IF (field%ndim/=2 .OR. field%data_type/=type_integer) STOP 'get_val_i2d : bad pointer assign ation'271 IF (field%ndim/=2 .OR. field%data_type/=type_integer) STOP 'get_val_i2d : bad pointer assignment' 271 272 field_pt=>field%ival2d 272 273 END SUBROUTINE getval_i2d … … 277 278 TYPE(t_field),INTENT(IN) :: field 278 279 279 IF (field%ndim/=3 .OR. field%data_type/=type_integer) STOP 'get_val_i3d : bad pointer assign ation'280 IF (field%ndim/=3 .OR. field%data_type/=type_integer) STOP 'get_val_i3d : bad pointer assignment' 280 281 field_pt=>field%ival3d 281 282 END SUBROUTINE getval_i3d … … 286 287 TYPE(t_field),INTENT(IN) :: field 287 288 288 IF (field%ndim/=4 .OR. field%data_type/=type_integer) STOP 'get_val_i4d : bad pointer assign ation'289 IF (field%ndim/=4 .OR. field%data_type/=type_integer) STOP 'get_val_i4d : bad pointer assignment' 289 290 field_pt=>field%ival4d 290 291 END SUBROUTINE getval_i4d … … 295 296 TYPE(t_field),INTENT(IN) :: field 296 297 297 IF (field%ndim/=2 .OR. field%data_type/=type_logical) STOP 'get_val_l2d : bad pointer assign ation'298 IF (field%ndim/=2 .OR. field%data_type/=type_logical) STOP 'get_val_l2d : bad pointer assignment' 298 299 field_pt=>field%lval2d 299 300 END SUBROUTINE getval_l2d … … 304 305 TYPE(t_field),INTENT(IN) :: field 305 306 306 IF (field%ndim/=3 .OR. field%data_type/=type_logical) STOP 'get_val_l3d : bad pointer assign ation'307 IF (field%ndim/=3 .OR. field%data_type/=type_logical) STOP 'get_val_l3d : bad pointer assignment' 307 308 field_pt=>field%lval3d 308 309 END SUBROUTINE getval_l3d … … 313 314 TYPE(t_field),INTENT(IN) :: field 314 315 315 IF (field%ndim/=4 .OR. field%data_type/=type_logical) STOP 'get_val_l4d : bad pointer assign ation'316 IF (field%ndim/=4 .OR. field%data_type/=type_logical) STOP 'get_val_l4d : bad pointer assignment' 316 317 field_pt=>field%lval4d 317 318 END SUBROUTINE getval_l4d -
codes/icosagcm/trunk/src/timeloop_gcm.f90
r158 r159 11 11 TYPE(t_message) :: req_ps0, req_theta_rhodz0, req_u0, req_q0 12 12 13 TYPE(t_field),POINTER :: f_phis(:)14 13 TYPE(t_field),POINTER :: f_q(:) 15 TYPE(t_field),POINTER :: f_dtheta(:), f_rhodz(:) 16 TYPE(t_field),POINTER :: f_ps(:),f_psm1(:), f_psm2(:) 17 TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:) 18 TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:) 19 TYPE(t_field),POINTER :: f_dps(:),f_dpsm1(:), f_dpsm2(:) 20 TYPE(t_field),POINTER :: f_du(:),f_dum1(:),f_dum2(:) 21 TYPE(t_field),POINTER :: f_dtheta_rhodz(:),f_dtheta_rhodzm1(:),f_dtheta_rhodzm2(:) 14 TYPE(t_field),POINTER :: f_rhodz(:), f_mass(:), f_massm1(:), f_massm2(:), f_dmass(:) 15 TYPE(t_field),POINTER :: f_phis(:), f_ps(:),f_psm1(:), f_psm2(:), f_dps(:) 16 TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:), f_du(:) 17 TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:), f_dtheta_rhodz(:) 22 18 TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:) 23 19 … … 34 30 USE caldyn_mod 35 31 USE etat0_mod 32 USE disvert_mod 36 33 USE guided_mod 37 34 USE advect_tracer_mod … … 45 42 IMPLICIT NONE 46 43 47 CHARACTER(len=255) :: scheme_name 48 49 CALL allocate_field(f_dps,field_t,type_real) 50 CALL allocate_field(f_du,field_u,type_real,llm) 51 CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm) 52 ! Model state at current time step (RK/MLF/Euler) 53 CALL allocate_field(f_phis,field_t,type_real) 54 CALL allocate_field(f_ps,field_t,type_real) 55 CALL allocate_field(f_u,field_u,type_real,llm) 56 CALL allocate_field(f_theta_rhodz,field_t,type_real,llm) 57 ! Model state at previous time step (RK/MLF) 58 CALL allocate_field(f_psm1,field_t,type_real) 59 CALL allocate_field(f_um1,field_u,type_real,llm) 60 CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm) 61 ! Tracers 62 CALL allocate_field(f_q,field_t,type_real,llm,nqtot) 63 CALL allocate_field(f_rhodz,field_t,type_real,llm) 64 ! Mass fluxes 65 CALL allocate_field(f_hflux,field_u,type_real,llm) ! instantaneous mass fluxes 66 CALL allocate_field(f_wflux,field_t,type_real,llm+1) ! computed by caldyn 67 CALL allocate_field(f_hfluxt,field_u,type_real,llm) ! mass "fluxes" accumulated in time 68 CALL allocate_field(f_wfluxt,field_t,type_real,llm+1) 69 44 CHARACTER(len=255) :: def 70 45 71 46 !---------------------------------------------------- … … 122 97 123 98 124 125 126 127 scheme_name='runge_kutta' 128 CALL getin('scheme',scheme_name) 129 130 SELECT CASE (TRIM(scheme_name)) 99 def='eta_mass' 100 CALL getin('caldyn_eta',def) 101 SELECT CASE(TRIM(def)) 102 CASE('eta_mass') 103 caldyn_eta=eta_mass 104 CASE('eta_lag') 105 caldyn_eta=eta_lag 106 CASE DEFAULT 107 IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_eta : <', TRIM(def),'> options are <eta_mass>, <eta_lag>' 108 STOP 109 END SELECT 110 IF (is_mpi_root) PRINT *, 'caldyn_eta=',TRIM(def), caldyn_eta, eta_lag, eta_mass 111 112 ! Time-independant orography 113 CALL allocate_field(f_phis,field_t,type_real,name='phis') 114 ! Trends 115 CALL allocate_field(f_du,field_u,type_real,llm,name='du') 116 CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm,name='dtheta_rhodz') 117 ! Model state at current time step (RK/MLF/Euler) 118 CALL allocate_field(f_ps,field_t,type_real, name='ps') 119 CALL allocate_field(f_mass,field_t,type_real,llm,name='mass') 120 CALL allocate_field(f_u,field_u,type_real,llm,name='u') 121 CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,name='theta_rhodz') 122 ! Model state at previous time step (RK/MLF) 123 CALL allocate_field(f_um1,field_u,type_real,llm,name='um1') 124 CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm,name='theta_rhodzm1') 125 ! Tracers 126 CALL allocate_field(f_q,field_t,type_real,llm,nqtot) 127 CALL allocate_field(f_rhodz,field_t,type_real,llm,name='rhodz') 128 ! Mass fluxes 129 CALL allocate_field(f_hflux,field_u,type_real,llm) ! instantaneous mass fluxes 130 CALL allocate_field(f_hfluxt,field_u,type_real,llm) ! mass "fluxes" accumulated in time 131 CALL allocate_field(f_wflux,field_t,type_real,llm+1) ! vertical mass fluxes 132 133 IF(caldyn_eta == eta_mass) THEN ! eta = mass coordinate (default) 134 CALL allocate_field(f_dps,field_t,type_real,name='dps') 135 CALL allocate_field(f_psm1,field_t,type_real,name='psm1') 136 CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt') 137 ! the following are unused but must point to something 138 f_massm1 => f_mass 139 f_dmass => f_mass 140 ELSE ! eta = Lagrangian vertical coordinate 141 CALL allocate_field(f_mass,field_t,type_real,llm) 142 CALL allocate_field(f_massm1,field_t,type_real,llm) 143 CALL allocate_field(f_dmass,field_t,type_real,llm) 144 ! the following are unused but must point to something 145 f_wfluxt => f_wflux 146 f_dps => f_phis 147 f_psm1 => f_phis 148 END IF 149 150 151 def='runge_kutta' 152 CALL getin('scheme',def) 153 154 SELECT CASE (TRIM(def)) 131 155 CASE('euler') 132 156 scheme=euler … … 141 165 nb_stage=matsuno_period+1 142 166 ! Model state 2 time steps ago (MLF) 143 CALL allocate_field(f_psm2,field_t,type_real)144 167 CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm) 145 168 CALL allocate_field(f_um2,field_u,type_real,llm) 169 IF(caldyn_eta == eta_mass) THEN ! eta = mass coordinate (default) 170 CALL allocate_field(f_psm2,field_t,type_real) 171 ! the following are unused but must point to something 172 f_massm2 => f_mass 173 ELSE ! eta = Lagrangian vertical coordinate 174 CALL allocate_field(f_massm2,field_t,type_real,llm) 175 ! the following are unused but must point to something 176 f_psm2 => f_phis 177 END IF 178 146 179 CASE default 147 PRINT*,'Bad selector for variable scheme : <', TRIM( scheme_name), &180 PRINT*,'Bad selector for variable scheme : <', TRIM(def), & 148 181 ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>' 149 182 STOP … … 157 190 CALL init_check_conserve 158 191 CALL init_physics 159 CALL etat0(f_ps,f_ phis,f_theta_rhodz,f_u, f_q)192 CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q) 160 193 161 194 CALL transfert_request(f_phis,req_i0) … … 173 206 USE icosa 174 207 USE dissip_gcm_mod 208 USE disvert_mod 175 209 USE caldyn_mod 210 USE caldyn_gcm_mod, ONLY : req_ps 176 211 USE etat0_mod 177 212 USE guided_mod … … 184 219 USE check_conserve_mod 185 220 IMPLICIT NONE 186 REAL(rstd),POINTER :: phis(:)187 221 REAL(rstd),POINTER :: q(:,:,:) 188 REAL(rstd),POINTER :: ps(:) ,psm1(:), psm2(:) 189 REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:) 190 REAL(rstd),POINTER :: rhodz(:,:), theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:) 191 REAL(rstd),POINTER :: dps(:), dpsm1(:), dpsm2(:) 192 REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:) 193 REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:) 222 REAL(rstd),POINTER :: phis(:), ps(:) ,psm1(:), psm2(:), dps(:) 223 REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:), du(:,:) 224 REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), massm1(:,:), massm2(:,:), dmass(:,:) 225 REAL(rstd),POINTER :: theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:), dtheta_rhodz(:,:) 194 226 REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:) 195 227 196 228 INTEGER :: ind 197 229 INTEGER :: it,i,j,n, stage 198 CHARACTER(len=255) :: scheme_name199 230 LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time 200 231 LOGICAL, PARAMETER :: check=.FALSE. 201 232 202 CALL caldyn_BC(f_phis, f_wflux) 233 CALL caldyn_BC(f_phis, f_wflux) ! set constant values in first/last interfaces 203 234 204 235 !$OMP BARRIER … … 234 265 DO stage=1,nb_stage 235 266 CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), & 236 f_phis,f_ps,f_ theta_rhodz,f_u, f_q, &267 f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q, & 237 268 f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) 238 269 SELECT CASE (scheme) … … 333 364 334 365 SUBROUTINE RK_scheme(stage) 335 USE caldyn_gcm_mod336 366 IMPLICIT NONE 337 367 INTEGER :: ind, stage … … 344 374 tau = dt*coef(stage) 345 375 346 DO ind=1,ndomain 347 CALL swap_dimensions(ind) 348 CALL swap_geometry(ind) 349 ps=f_ps(ind) 350 psm1=f_psm1(ind) 351 dps=f_dps(ind) 352 353 IF (stage==1) THEN ! first stage : save model state in XXm1 354 IF (omp_first) THEN 355 DO j=jj_begin,jj_end 356 DO i=ii_begin,ii_end 357 ij=(j-1)*iim+i 358 psm1(ij)=ps(ij) 376 ! if mass coordinate, deal with ps first on one core 377 IF(caldyn_eta==eta_mass) THEN 378 IF (omp_first) THEN 379 DO ind=1,ndomain 380 CALL swap_dimensions(ind) 381 CALL swap_geometry(ind) 382 ps=f_ps(ind) 383 psm1=f_psm1(ind) 384 dps=f_dps(ind) 385 386 IF (stage==1) THEN ! first stage : save model state in XXm1 387 DO j=jj_begin,jj_end 388 DO i=ii_begin,ii_end 389 ij=(j-1)*iim+i 390 psm1(ij)=ps(ij) 391 ENDDO 392 ENDDO 393 ENDIF 394 395 ! updates are of the form x1 := x0 + tau*f(x1) 396 DO j=jj_begin,jj_end 397 DO i=ii_begin,ii_end 398 ij=(j-1)*iim+i 399 ps(ij)=psm1(ij)+tau*dps(ij) 400 ENDDO 359 401 ENDDO 360 ENDDO 361 ENDIF 362 END IF 363 364 ! updates are of the form x1 := x0 + tau*f(x1) 365 IF (omp_first) THEN 366 DO j=jj_begin,jj_end 367 DO i=ii_begin,ii_end 368 ij=(j-1)*iim+i 369 ps(ij)=psm1(ij)+tau*dps(ij) 370 ENDDO 371 ENDDO 372 ENDIF 373 374 ENDDO 375 376 CALL send_message(f_ps,req_ps) 377 378 379 402 ENDDO 403 ENDIF 404 405 CALL send_message(f_ps,req_ps) 406 END IF 407 408 ! now deal with other prognostic variables 380 409 DO ind=1,ndomain 381 410 CALL swap_dimensions(ind) … … 386 415 387 416 IF (stage==1) THEN ! first stage : save model state in XXm1 388 389 DO l=ll_begin,ll_end417 418 DO l=ll_begin,ll_end 390 419 DO j=jj_begin,jj_end 391 420 DO i=ii_begin,ii_end … … 398 427 ENDDO 399 428 ENDDO 429 430 IF(caldyn_eta==eta_lag) THEN ! mass = additional prognostic variable 431 DO l=ll_begin,ll_end 432 DO j=jj_begin,jj_end 433 DO i=ii_begin,ii_end 434 ij=(j-1)*iim+i 435 massm1(ij,l)=mass(ij,l) 436 ENDDO 437 ENDDO 438 ENDDO 439 END IF 400 440 401 441 END IF … … 413 453 ENDDO 414 454 ENDDO 415 455 IF(caldyn_eta==eta_lag) THEN ! mass = additional prognostic variable 456 DO l=ll_begin,ll_end 457 DO j=jj_begin,jj_end 458 DO i=ii_begin,ii_end 459 ij=(j-1)*iim+i 460 mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l) 461 ENDDO 462 ENDDO 463 ENDDO 464 END IF 465 416 466 IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage 417 467 hflux=f_hflux(ind); hfluxt=f_hfluxt(ind) … … 476 526 END SUBROUTINE timeloop 477 527 478 SUBROUTINE compute_rhodz(comp, ps, rhodz) 479 USE icosa 480 USE disvert_mod 481 USE omp_para 482 LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check 483 REAL(rstd), INTENT(IN) :: ps(iim*jjm) 484 REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm) 485 REAL(rstd) :: m, err 486 INTEGER :: l,i,j,ij,dd 487 err=0. 488 489 dd=0 490 491 DO l = ll_begin, ll_end 492 DO j=jj_begin-dd,jj_end+dd 493 DO i=ii_begin-dd,ii_end+dd 494 ij=(j-1)*iim+i 495 m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 496 IF(comp) THEN 497 rhodz(ij,l) = m 498 ELSE 499 err = MAX(err,abs(m-rhodz(ij,l))) 500 END IF 501 ENDDO 502 ENDDO 503 ENDDO 504 505 IF(.NOT. comp) THEN 506 IF(err>1e-10) THEN 507 PRINT *, 'Discrepancy between ps and rhodz detected', err 508 STOP 509 ELSE 510 ! PRINT *, 'No discrepancy between ps and rhodz detected' 511 END IF 512 END IF 513 514 END SUBROUTINE compute_rhodz 515 516 SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero) 528 SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero) 517 529 USE icosa 518 530 USE omp_para 531 USE disvert_mod 532 IMPLICIT NONE 519 533 REAL(rstd), INTENT(IN) :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1) 520 534 REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1) … … 538 552 ENDDO 539 553 540 DO l=ll_begin,ll_endp1 541 DO j=jj_begin,jj_end 542 DO i=ii_begin,ii_end 543 ij=(j-1)*iim+i 544 wfluxt(ij,l) = tau*wflux(ij,l) 545 ENDDO 546 ENDDO 547 ENDDO 554 IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag 555 DO l=ll_begin,ll_endp1 556 DO j=jj_begin,jj_end 557 DO i=ii_begin,ii_end 558 ij=(j-1)*iim+i 559 wfluxt(ij,l) = tau*wflux(ij,l) 560 ENDDO 561 ENDDO 562 ENDDO 563 END IF 548 564 ELSE 549 565 … … 559 575 ENDDO 560 576 561 DO l=ll_begin,ll_endp1 562 DO j=jj_begin,jj_end 563 DO i=ii_begin,ii_end 564 ij=(j-1)*iim+i 565 wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l) 566 ENDDO 567 ENDDO 568 ENDDO 569 570 END IF 577 IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag 578 DO l=ll_begin,ll_endp1 579 DO j=jj_begin,jj_end 580 DO i=ii_begin,ii_end 581 ij=(j-1)*iim+i 582 wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l) 583 ENDDO 584 ENDDO 585 ENDDO 586 END IF 587 588 END IF 571 589 572 590 END SUBROUTINE accumulate_fluxes
Note: See TracChangeset
for help on using the changeset viewer.