![(please configure the [header_logo] section in trac.ini)](/dynamico/chrome/site/your_project_logo.png){kind=logo}
Changeset 196 for codes/icosagcm/trunk/src/physics_dcmip.f90
- Timestamp:
- 07/06/14 09:10:13 (10 years ago)
- File:
-
- 1 edited
-
codes/icosagcm/trunk/src/physics_dcmip.f90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/physics_dcmip.f90
r192 r196 6 6 !$OMP THREADPRIVATE(testcase) 7 7 8 PUBLIC init_physics, physics9 8 TYPE(t_field),POINTER :: f_out_i(:) 10 9 TYPE(t_field),POINTER :: f_precl(:) 11 REAL(rstd),POINTER :: out_i(:,:) 10 REAL(rstd),POINTER :: out_i(:,:) 11 12 PUBLIC :: compute_phys_wrap, init_physics 12 13 13 14 CONTAINS 14 15 15 16 SUBROUTINE init_physics 16 IMPLICIT NONE 17 17 IMPLICIT NONE 18 18 testcase=1 ! OK for 4.2 (moist baroclinic instability) 19 19 CALL getin("dcmip_physics",testcase) 20 CALL allocate_field(f_out_i,field_t,type_real,llm)21 CALL allocate_field(f_precl,field_t,type_real)22 23 20 END SUBROUTINE init_physics 24 21 25 26 SUBROUTINE physics( it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q) 27 USE icosa 28 IMPLICIT NONE 29 INTEGER,INTENT(IN) :: it 30 TYPE(t_field),POINTER :: f_phis(:) 31 TYPE(t_field),POINTER :: f_ps(:) 32 TYPE(t_field),POINTER :: f_theta_rhodz(:) 33 TYPE(t_field),POINTER :: f_ue(:) 34 TYPE(t_field),POINTER :: f_q(:) 35 36 REAL(rstd),POINTER :: phis(:) 37 REAL(rstd),POINTER :: ps(:) 38 REAL(rstd),POINTER :: theta_rhodz(:,:) 39 REAL(rstd),POINTER :: ue(:,:) 40 REAL(rstd),POINTER :: q(:,:,:) 41 REAL(rstd),POINTER :: precl(:) 42 INTEGER :: ind 43 44 CALL transfert_request(f_ue,req_e1_vect) 45 DO ind=1,ndomain 46 IF (.NOT. assigned_domain(ind)) CYCLE 47 CALL swap_dimensions(ind) 48 CALL swap_geometry(ind) 49 50 phis=f_phis(ind) 51 ps=f_ps(ind) 52 theta_rhodz=f_theta_rhodz(ind) 53 ue=f_ue(ind) 54 q=f_q(ind) 55 out_i=f_out_i(ind) 56 precl=f_precl(ind) 57 58 CALL compute_physics( phis, ps, theta_rhodz, ue, q(:,:,1), precl) 59 60 ENDDO 61 62 ! CALL writefield("out_i",f_out_i) 63 64 IF (mod(it,itau_out)==0 ) THEN 65 CALL writefield("precl",f_precl) 66 ENDIF 67 68 END SUBROUTINE physics 69 70 SUBROUTINE compute_physics(phis, ps, theta_rhodz, ue, q, precl) 71 USE icosa 72 USE pression_mod 73 USE exner_mod 74 USE theta2theta_rhodz_mod 75 USE geopotential_mod 76 USE wind_mod 77 IMPLICIT NONE 78 REAL(rstd) :: phis(iim*jjm) 79 REAL(rstd) :: ps(iim*jjm) 80 REAL(rstd) :: theta_rhodz(iim*jjm,llm) 81 REAL(rstd) :: ue(3*iim*jjm,llm) 82 REAL(rstd) :: q(iim*jjm,llm) 83 REAL(rstd) :: precl(iim*jjm) 84 85 REAL(rstd) :: p(iim*jjm,llm+1) 86 REAL(rstd) :: pks(iim*jjm) 87 REAL(rstd) :: pk(iim*jjm,llm) 88 REAL(rstd) :: phi(iim*jjm,llm) 89 REAL(rstd) :: T(iim*jjm,llm) 90 REAL(rstd) :: Tfi(iim*jjm,llm) 91 REAL(rstd) :: theta(iim*jjm,llm) 92 93 REAL(rstd) :: uc(iim*jjm,3,llm) 94 REAL(rstd) :: u(iim*jjm,llm) 95 REAL(rstd) :: v(iim*jjm,llm) 96 REAL(rstd) :: ufi(iim*jjm,llm) 97 REAL(rstd) :: vfi(iim*jjm,llm) 98 REAL(rstd) :: qfi(iim*jjm,llm) 99 REAL(rstd) :: utemp(iim*jjm,llm) 100 REAL(rstd) :: vtemp(iim*jjm,llm) 101 REAL(rstd) :: lat(iim*jjm) 102 REAL(rstd) :: lon 103 REAL(rstd) :: pmid(iim*jjm,llm) 104 REAL(rstd) :: pint(iim*jjm,llm+1) 105 REAL(rstd) :: pdel(iim*jjm,llm) 106 INTEGER :: i,j,l,ij 107 108 ! PRINT *,'Entering in DCMIP physics' 109 CALL compute_pression(ps,p,0) 110 CALL compute_exner(ps,p,pks,pk,0) 111 CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,0) 112 CALL compute_geopotential(phis,pks,pk,theta,phi,0) 113 CALL compute_theta_rhodz2temperature(ps,theta_rhodz,T,0) 114 ! Reconstruct wind vector at hexagons 115 CALL compute_wind_centered(ue,uc) 116 CALL compute_wind_centered_lonlat_compound(uc, u, v) 117 118 ! Convert from Tv to T 119 DO l=1,llm 120 DO j=jj_begin,jj_end 121 DO i=ii_begin,ii_end 122 ij=(j-1)*iim+i 123 T(ij,l)=T(ij,l)/(1+0.608*q(ij,l)) 124 ENDDO 125 ENDDO 126 ENDDO 127 128 DO j=jj_begin,jj_end 129 DO i=ii_begin,ii_end 130 ij=(j-1)*iim+i 131 CALL xyz2lonlat(xyz_i(ij,:),lon,lat(ij)) 132 ENDDO 133 ENDDO 134 135 ! bottom-up indexing (DYNAMICO) : u,utemp, v,vtemp 136 ! top-down vertical indexing (DCMIP) : ufi, vfi, ... 137 ! => copy fields and mirror vertical index 22 SUBROUTINE compute_phys_wrap(args) 23 USE physics_interface_mod 24 TYPE(physics_inout) :: args 25 CALL compute_physics(args%ngrid, args%dt_phys, args%lat, & 26 args%p, args%Temp, args%ulon, args%ulat, args%q(:,:,1), & 27 args%dTemp, args%dulon, args%dulat, args%dq(:,:,1), args%extra_2D(:,1)) 28 END SUBROUTINE compute_phys_wrap 29 30 SUBROUTINE compute_physics(ngrid,dt_phys,lat, p,Temp,u,v,q, dTemp,du,dv,dq, precl) 31 USE icosa 32 IMPLICIT NONE 33 INTEGER :: ngrid 34 REAL(rstd) :: lat(ngrid) 35 REAL(rstd) :: ps(ngrid) 36 REAL(rstd) :: precl(ngrid) 37 ! arguments with bottom-up indexing (DYNAMICO) 38 REAL(rstd) :: p(ngrid,llm+1) 39 REAL(rstd) :: Temp(ngrid,llm) 40 REAL(rstd) :: u(ngrid,llm) 41 REAL(rstd) :: v(ngrid,llm) 42 REAL(rstd) :: q(ngrid,llm) 43 REAL(rstd) :: dTemp(ngrid,llm) 44 REAL(rstd) :: du(ngrid,llm) 45 REAL(rstd) :: dv(ngrid,llm) 46 REAL(rstd) :: dq(ngrid,llm) 47 ! local arrays with top-down vertical indexing (DCMIP) 48 REAL(rstd) :: pint(ngrid,llm+1) 49 REAL(rstd) :: pmid(ngrid,llm) 50 REAL(rstd) :: pdel(ngrid,llm) 51 REAL(rstd) :: Tfi(ngrid,llm) 52 REAL(rstd) :: ufi(ngrid,llm) 53 REAL(rstd) :: vfi(ngrid,llm) 54 REAL(rstd) :: qfi(ngrid,llm) 55 56 INTEGER :: i,j,l,ll,ij 57 REAL(rstd) :: dt_phys, inv_dt 58 59 ! prepare input fields and mirror vertical index 60 ps(:) = p(:,1) ! surface pressure 61 138 62 DO l=1,llm+1 139 63 DO j=jj_begin,jj_end … … 144 68 ENDDO 145 69 ENDDO 146 147 ! Pressure inside layers 70 148 71 DO l=1,llm 149 DO j=jj_begin,jj_end 150 DO i=ii_begin,ii_end 151 ij=(j-1)*iim+i 152 pmid(ij,l)=0.5*(pint(ij,l)+pint(ij,l+1)) 153 ENDDO 154 ENDDO 72 ll=llm+1-l 73 DO j=jj_begin,jj_end 74 DO i=ii_begin,ii_end 75 ij=(j-1)*iim+i 76 pmid(ij,l)=0.5*(pint(ij,l)+pint(ij,l+1)) ! Pressure inside layers 77 pdel(ij,l)=pint(ij,l+1)-pint(ij,l) ! Pressure difference between two layers 78 ufi(ij,l)=u(ij,ll) 79 vfi(ij,l)=v(ij,ll) 80 qfi(ij,l)=q(ij,ll) 81 Tfi(ij,l)=Temp(ij,ll)/(1+0.608*qfi(ij,l)) 82 ENDDO 83 ENDDO 155 84 ENDDO 156 85 157 ! Pressure difference between two layers 86 CALL simple_physics(ngrid, llm, dt_phys, lat, tfi, qfi , ufi, vfi, pmid, pint, pdel, 1./pdel, ps, precl, testcase) 87 88 ! Mirror vertical index and compute tendencies 89 inv_dt = 1./dt_phys 158 90 DO l=1,llm 159 DO j=jj_begin,jj_end 160 DO i=ii_begin,ii_end 161 ij=(j-1)*iim+i 162 pdel(ij,l)=pint(ij,l+1)-pint(ij,l) 163 ENDDO 164 ENDDO 165 ENDDO 166 167 ! Copy T,u,v,q for input to physics 168 DO l=1,llm 169 DO j=jj_begin,jj_end 170 DO i=ii_begin,ii_end 171 ij=(j-1)*iim+i 172 Tfi(ij,l)=T(ij,llm+1-l) 173 ufi(ij,l)=u(ij,llm+1-l) 174 vfi(ij,l)=v(ij,llm+1-l) 175 qfi(ij,l)=q(ij,llm+1-l) 176 ENDDO 177 ENDDO 178 ENDDO 179 180 CALL simple_physics(iim*jjm, llm, dt, lat, tfi, qfi , ufi, vfi, pmid, pint, pdel, 1/pdel, ps, precl, testcase) 181 182 ! Copy back T,u,v,q and mirror vertical index 183 DO l=1,llm 184 DO j=jj_begin,jj_end 185 DO i=ii_begin,ii_end 186 ij=(j-1)*iim+i 187 T(ij,l)=Tfi(ij,llm+1-l) 188 utemp(ij,l)=ufi(ij,llm+1-l) 189 vtemp(ij,l)=vfi(ij,llm+1-l) 190 q(ij,l)=qfi(ij,llm+1-l) 191 ENDDO 192 ENDDO 193 ENDDO 194 195 196 ! Convert back T to Tv 197 DO l=1,llm 198 DO j=jj_begin,jj_end 199 DO i=ii_begin,ii_end 200 ij=(j-1)*iim+i 201 T(ij,l)=T(ij,l)*(1+0.608*q(ij,l)) 202 ENDDO 203 ENDDO 204 ENDDO 205 206 ! Compute velocity update at hexagons 207 utemp=utemp-u 208 vtemp=vtemp-v 209 210 ! lon-lat -> 3D 211 DO l=1,llm 212 DO j=jj_begin,jj_end 213 DO i=ii_begin,ii_end 214 ij=(j-1)*iim+i 215 uc(ij,:,l)=utemp(ij,l)*elon_i(ij,:)+vtemp(ij,l)*elat_i(ij,:) 216 ENDDO 217 ENDDO 91 ll=llm+1-l 92 DO j=jj_begin,jj_end 93 DO i=ii_begin,ii_end 94 ij=(j-1)*iim+i 95 dTemp(ij,l) = inv_dt * ( Tfi(ij,ll)*(1+0.608*qfi(ij,ll)) - Temp(ij,l) ) 96 du(ij,l) = inv_dt * (ufi(ij,ll) - u(ij,l)) 97 dv(ij,l) = inv_dt * (vfi(ij,ll) - v(ij,l)) 98 dq(ij,l) = inv_dt * (qfi(ij,ll) - q(ij,l)) 99 ENDDO 100 ENDDO 218 101 ENDDO 219 220 ! Update velocity at velocity points221 DO l=1,llm222 DO j=jj_begin,jj_end223 DO i=ii_begin,ii_end224 ij=(j-1)*iim+i225 ue(ij+u_right,l)=ue(ij+u_right,l)+sum( 0.5*(uc(ij,:,l) + uc(ij+t_right,:,l))*ep_e(ij+u_right,:) )226 ue(ij+u_lup,l)=ue(ij+u_lup,l)+sum( 0.5*(uc(ij,:,l) + uc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )227 ue(ij+u_ldown,l)=ue(ij+u_ldown,l)+sum( 0.5*(uc(ij,:,l) + uc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )228 ENDDO229 ENDDO230 ENDDO231 232 CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0)233 234 102 235 103 END SUBROUTINE compute_physics 236 104 105 END MODULE physics_dcmip_mod 237 106 238 107 subroutine simple_physics (pcols, pver, dtime, lat, t, q, u, v, pmid, pint, pdel, rpdel, ps, precl, test) … … 501 370 do i=1,pcols 502 371 qsat = epsilo*e0/pmid(i,k)*exp(-latvap/rh2o*((1._r8/t(i,k))-1._r8/T0)) ! saturation specific humidity 503 out_i(i,llm+1-k)=q(i,k)-qsat372 ! out_i(i,llm+1-k)=q(i,k)-qsat 504 373 if (q(i,k) > qsat) then ! saturated? 505 374 tmp = 1._r8/dtime*(q(i,k)-qsat)/(1._r8+(latvap/cpair)*(epsilo*latvap*qsat/(rair*t(i,k)**2))) … … 694 563 end subroutine simple_physics 695 564 696 697 698 699 END MODULE physics_dcmip_mod
Note: See TracChangeset
for help on using the changeset viewer.
