Changeset 48 for codes/icosagcm/trunk
- Timestamp:
- 07/29/12 03:14:26 (12 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 1 added
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/etat0_dcmip3.f90
r37 r48 7 7 8 8 USE genmod 9 USE dcmip_initial_conditions_test_1_2_3 10 9 11 PRIVATE 10 12 ! some of the following should be obtained via getin … … 15 17 REAL(rstd),PARAMETER :: u0=20. ! Maximum amplitude of the zonal wind (m.s-1) 16 18 REAL(rstd),PARAMETER :: zs=0. ! Surface elevation (m) 17 REAL(rstd),PARAMETER :: N=0.01 ! Brunt-Vais sala frequency (s-1)19 REAL(rstd),PARAMETER :: N=0.01 ! Brunt-Vaisala frequency (s-1) 18 20 REAL(rstd),PARAMETER :: Teq=300. ! Surface temperature at the equator (K) 19 21 REAL(rstd),PARAMETER :: Peq=1e5 ! Reference surface pressure at the equator (hPa) … … 74 76 REAL(rstd) :: Ts(iim*jjm) 75 77 REAL(rstd) :: s(iim*jjm) 76 REAL(rstd) :: z(iim*jjm,llm)77 78 REAL(rstd) :: T(iim*jjm,llm) 78 79 REAL(rstd) :: p(iim*jjm,llm+1) 79 REAL(rstd) :: thetap(iim*jjm,llm) 80 REAL(rstd) :: thetap_rhodz(iim*jjm,llm) 80 REAL(rstd) :: theta(iim*jjm,llm) 81 81 REAL(rstd) :: ulon(3*iim*jjm,llm) 82 82 REAL(rstd) :: ulat(3*iim*jjm,llm) … … 86 86 REAL(rstd) :: Rd ! gas constant of dry air, P=rho.Rd.T 87 87 REAL(rstd) :: alpha, rm, zs 88 REAL(rstd) :: lon,lat 89 REAL(rstd) :: C0,C1 90 REAL(rstd) :: GG 91 REAL(rstd) :: pspsk,r 92 93 Rd=cpp/kappa 94 95 ! alpha=g/Rd/lapse_rate ! g/(Rd*Gamma) 96 97 GG=(g/N)**2/cpp 98 C0=0.25*u0*(u0+2.*Omega*radius) 99 100 q(:,:,:)=0 101 102 103 DO j=jj_begin,jj_end 104 DO i=ii_begin,ii_end 88 REAL(rstd) :: lon,lat, C0, C1, GG 89 REAL(rstd) :: p0psk, pspsk,r,zz, thetab, thetap 90 REAL(rstd) :: dummy, pp 91 LOGICAL :: use_dcmip_routine 92 93 Rd=cpp*kappa 94 95 GG=(g/N)**2/cpp 96 C0=0.25*u0*(u0+2.*Omega*radius) 97 98 q(:,:,:)=0 99 100 ! use_dcmip_routine=.TRUE. 101 use_dcmip_routine=.FALSE. 102 dummy=0. 103 104 pp=peq 105 DO j=jj_begin,jj_end 106 DO i=ii_begin,ii_end 105 107 ij=(j-1)*iim+i 106 108 CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 107 C1=C0*(cos(2*lat)-1) 108 109 !--- GROUND GEOPOTENTIAL 110 phis(ij)=0. 111 112 !--- GROUND TEMPERATURE 113 Ts(ij) = GG+(Teq-GG)*EXP(-C1*(N/g)**2) 114 115 !--- GROUND PRESSURE 116 Ps(ij) = peq*EXP(C1/GG/Rd)*(Ts(ij)/Teq)**(1/kappa) 117 118 119 r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) 120 s(ij)= d**2/(d**2+r**2) 121 END DO 122 END DO 123 124 CALL compute_pression(ps,p,0) 125 126 DO l=1,llm 127 DO j=jj_begin,jj_end 109 110 IF(use_dcmip_routine) THEN 111 CALL test3_gravity_wave(lon,lat,pp,dummy,0, dummy,dummy,dummy,dummy,phis(ij),ps(ij),dummy,dummy) 112 ELSE 113 C1=C0*(cos(2*lat)-1) 114 115 !--- GROUND GEOPOTENTIAL 116 phis(ij)=0. 117 118 !--- GROUND TEMPERATURE 119 Ts(ij) = GG+(Teq-GG)*EXP(-C1*(N/g)**2) 120 121 !--- GROUND PRESSURE 122 Ps(ij) = peq*EXP(C1/GG/Rd)*(Ts(ij)/Teq)**(1/kappa) 123 124 125 r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) 126 s(ij)= d**2/(d**2+r**2) 127 END IF 128 END DO 129 END DO 130 131 CALL compute_pression(ps,p,0) 132 133 DO l=1,llm 134 DO j=jj_begin,jj_end 128 135 DO i=ii_begin,ii_end 129 ij=(j-1)*iim+i 130 131 pspsk=(0.5*(p(ij,l+1)+p(ij,l))/ps(ij) )**kappa 132 T(ij,l)=Ts(ij)*pspsk / ( Ts(ij) / GG * ( pspsk-1) +1) 133 z(ij,l)=-g/N**2*log( Ts(ij)/GG * (pspsk -1)+1) 136 ij=(j-1)*iim+i 137 pp=0.5*(p(ij,l+1)+p(ij,l)) ! full-layer pressure 138 IF(use_dcmip_routine) THEN 139 CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 140 CALL test3_gravity_wave(lon,lat,pp,dummy,0,dummy,dummy,dummy,T(ij,l),dummy,dummy,dummy,dummy) 141 ELSE 142 pspsk=(pp/ps(ij))**kappa 143 p0psk=(Peq/ps(ij))**kappa 144 thetab = Ts(ij)*p0psk / ( Ts(ij) / GG * ( pspsk-1) +1) ! background pot. temp. 145 zz = -g/N**2*log( Ts(ij)/GG * (pspsk -1)+1) ! altitude 146 thetap = dtheta *sin(2*Pi*zz/Lz) * s(ij) ! perturbation pot. temp. 147 theta(ij,l) = thetab + thetap 148 T(ij,l) = theta(ij,l)* ((pp/peq)**kappa) 149 ! T(ij,l) = Ts(ij)*pspsk / ( Ts(ij) / GG * ( pspsk-1) +1) ! background temp. 150 END IF 134 151 ENDDO 135 ENDDO 136 ENDDO 137 138 CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0) 139 140 DO l=1,llm 141 DO j=jj_begin,jj_end 142 DO i=ii_begin,ii_end 143 ij=(j-1)*iim+i 144 thetap(ij,l)=dtheta *sin(2*Pi*z(ij,l)/Lz) * s(ij) 152 ENDDO 153 ENDDO 154 155 IF(use_dcmip_routine) THEN 156 CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0) 157 ELSE 158 ! CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0) 159 CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0) 160 END IF 161 162 pp=peq 163 DO l=1,llm 164 DO j=jj_begin-1,jj_end+1 165 DO i=ii_begin-1,ii_end+1 166 ij=(j-1)*iim+i 167 IF(use_dcmip_routine) THEN 168 CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat) 169 CALL test3_gravity_wave(lon,lat,pp,0.,0, & 170 ulon(ij+u_right,l),ulat(ij+u_right,l),& 171 dummy,dummy,dummy,dummy,dummy,dummy) 172 CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat) 173 CALL test3_gravity_wave(lon,lat,pp,0.,0,& 174 ulon(ij+u_lup,l),ulat(ij+u_lup,l),& 175 dummy,dummy,dummy,dummy,dummy,dummy) 176 CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat) 177 CALL test3_gravity_wave(lon,lat,pp,0.,0,& 178 ulon(ij+u_ldown,l),ulat(ij+u_ldown,l),& 179 dummy,dummy,dummy,dummy,dummy,dummy) 180 ELSE 181 CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat) 182 ulon(ij+u_right,l)=u0*cos(lat) 183 ulat(ij+u_right,l)=0 184 185 CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat) 186 ulon(ij+u_lup,l)=u0*cos(lat) 187 ulat(ij+u_lup,l)=0 188 189 CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat) 190 ulon(ij+u_ldown,l)=u0*cos(lat) 191 ulat(ij+u_ldown,l)=0 192 END IF 145 193 ENDDO 146 ENDDO 147 ENDDO 148 149 CALL compute_theta2theta_rhodz(ps,thetap,thetap_rhodz,0) 150 151 152 DO l=1,llm 153 DO j=jj_begin,jj_end 154 DO i=ii_begin,ii_end 155 ij=(j-1)*iim+i 156 theta_rhodz(ij,l)=theta_rhodz(ij,l)+thetap_rhodz(ij,l) 157 ENDDO 158 ENDDO 159 ENDDO 160 161 DO l=1,llm 162 DO j=jj_begin-1,jj_end+1 163 DO i=ii_begin-1,ii_end+1 164 ij=(j-1)*iim+i 165 CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat) 166 ulon(ij+u_right,l)=u0*cos(lat) 167 ulat(ij+u_right,l)=0 168 169 CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat) 170 ulon(ij+u_lup,l)=u0*cos(lat) 171 ulat(ij+u_lup,l)=0 172 173 CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat) 174 ulon(ij+u_ldown,l)=u0*cos(lat) 175 ulat(ij+u_ldown,l)=0 176 ENDDO 177 ENDDO 178 ENDDO 179 180 181 CALL compute_wind_perp_from_lonlat_compound(ulon,ulat,u) 182 183 184 185 END SUBROUTINE compute_etat0_DCMIP3 186 187 188 END MODULE etat0_DCMIP3_mod 194 ENDDO 195 ENDDO 196 197 CALL compute_wind_perp_from_lonlat_compound(ulon,ulat,u) 198 199 END SUBROUTINE compute_etat0_DCMIP3 200 201 202 END MODULE etat0_DCMIP3_mod
Note: See TracChangeset
for help on using the changeset viewer.