![(please configure the [header_logo] section in trac.ini)](/dynamico/chrome/site/your_project_logo.png){kind=logo}
Changeset 203 for codes/icosagcm/trunk
- Timestamp:
- 07/09/14 00:58:30 (10 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 4 edited
-
etat0.f90 (modified) (7 diffs)
-
etat0_dcmip5.f90 (modified) (4 diffs)
-
etat0_jablonowsky06.f90 (modified) (3 diffs)
-
vector.f90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/etat0.f90
r201 r203 15 15 USE mpipara, ONLY : is_mpi_root 16 16 USE disvert_mod 17 ! New interface 18 USE etat0_dcmip5_mod, ONLY : getin_etat0_dcmip5=>getin_etat0 19 ! Old interface 17 20 USE etat0_williamson_mod, ONLY : etat0_williamson_new 18 USE etat0_jablonowsky06_mod, ONLY : etat0_jablonowsky06=>etat019 21 USE etat0_academic_mod, ONLY : etat0_academic=>etat0 20 22 USE etat0_dcmip1_mod, ONLY : etat0_dcmip1=>etat0 … … 22 24 USE etat0_dcmip3_mod, ONLY : etat0_dcmip3=>etat0 23 25 USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>etat0 24 USE etat0_dcmip5_mod, ONLY : etat0_dcmip5=>etat025 26 USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0 26 27 USE dynetat0_gcm_mod, ONLY : dynetat0_start=>etat0 … … 51 52 52 53 SELECT CASE (TRIM(etat0_type)) 54 !------------------- New interface --------------------- 53 55 CASE ('isothermal') 54 56 CALL getin_etat0_isothermal 55 57 CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q) 58 CASE ('jablonowsky06') 59 CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q) 60 CASE ('dcmip5') 61 CALL getin_etat0_dcmip5 62 CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q) 63 64 !------------------- Old interface -------------------- 56 65 CASE ('williamson91.6') 57 66 init_mass=.FALSE. 58 67 CALL etat0_williamson_new(f_phis,f_mass,f_theta_rhodz,f_u, f_q) 59 CASE ('jablonowsky06')60 ! CALL etat0_jablonowsky06(f_ps,f_phis,f_theta_rhodz,f_u, f_q)61 CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)62 68 CASE ('academic') 63 69 CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q) … … 79 85 END IF 80 86 CALL etat0_dcmip4(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 81 CASE ('dcmip5')82 CALL etat0_dcmip5(f_ps,f_phis,f_theta_rhodz,f_u, f_q)83 87 CASE ('readnf_start') 84 88 print*,"readnf_start used" … … 135 139 q=f_q(ind) 136 140 CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, q) 137 138 141 ENDDO 139 142 END SUBROUTINE etat0_collocated … … 143 146 USE theta2theta_rhodz_mod 144 147 USE wind_mod 145 USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0_new 148 USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0 149 USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0 146 150 IMPLICIT NONE 147 151 REAL(rstd),INTENT(INOUT) :: ps(iim*jjm) … … 189 193 CALL compute_jablonowsky06(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i) 190 194 CALL compute_jablonowsky06(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e) 195 CASE('dcmip5') 196 CALL compute_dcmip5(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) 197 CALL compute_dcmip5(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) 191 198 END SELECT 192 199 -
codes/icosagcm/trunk/src/etat0_dcmip5.f90
r192 r203 16 16 REAL(rstd),PARAMETER :: epsilon=1e-25 17 17 REAL(rstd),PARAMETER :: Rd=287 18 REAL(rstd),SAVE :: lonc 19 !$OMP THREADPRIVATE(lonc) 20 REAL(rstd),SAVE :: latc 21 !$OMP THREADPRIVATE(latc) 22 TYPE(t_field),POINTER :: f_out_i(:) 23 REAL(rstd),POINTER,SAVE :: out_i(:,:) 24 !$OMP THREADPRIVATE(out_i) 18 19 REAL(rstd), PARAMETER :: lonc=Pi, latc=Pi/18, & 20 Tv0=T0*(1+0.608*q0), Tvt=Tv0-Gamma*zt 25 21 26 PUBLIC etat0 22 INTEGER, SAVE :: dcmip5_testcase 23 24 PUBLIC getin_etat0, compute_etat0 25 27 26 CONTAINS 28 27 29 30 31 SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 32 USE icosa 33 IMPLICIT NONE 34 TYPE(t_field),POINTER :: f_ps(:) 35 TYPE(t_field),POINTER :: f_phis(:) 36 TYPE(t_field),POINTER :: f_theta_rhodz(:) 37 TYPE(t_field),POINTER :: f_u(:) 38 TYPE(t_field),POINTER :: f_q(:) 39 40 REAL(rstd),POINTER :: ps(:) 41 REAL(rstd),POINTER :: phis(:) 42 REAL(rstd),POINTER :: theta_rhodz(:,:) 43 REAL(rstd),POINTER :: u(:,:) 44 REAL(rstd),POINTER :: q(:,:,:) 45 INTEGER :: ind 46 47 CALL allocate_field(f_out_i,field_t,type_real,llm) 48 49 DO ind=1,ndomain 50 IF (.NOT. assigned_domain(ind)) CYCLE 51 CALL swap_dimensions(ind) 52 CALL swap_geometry(ind) 53 ps=f_ps(ind) 54 phis=f_phis(ind) 55 theta_rhodz=f_theta_rhodz(ind) 56 u=f_u(ind) 57 q=f_q(ind) 58 out_i=f_out_i(ind) 59 CALL compute_etat0_dcmip5(ps, phis, theta_rhodz, u, q) 60 ENDDO 61 CALL writefield("out_i",f_out_i) 62 CALL deallocate_field(f_out_i) 63 END SUBROUTINE etat0 64 65 SUBROUTINE compute_etat0_dcmip5(ps, phis, theta_rhodz, ue, q) 66 USE icosa 67 USE pression_mod 68 USE wind_mod 69 USE theta2theta_rhodz_mod 70 USE exner_mod 71 IMPLICIT NONE 72 REAL(rstd),INTENT(OUT) :: ps(iim*jjm) 73 REAL(rstd),INTENT(OUT) :: phis(iim*jjm) 74 REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) 75 REAL(rstd),INTENT(OUT) :: ue(3*iim*jjm,llm) 76 REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot) 77 78 INTEGER :: i,j,l,ij 79 INTEGER :: testcase 80 REAL(rstd) :: Tv0, Tvt 81 REAL(rstd) :: r 82 REAL(rstd) :: z(iim*jjm,llm),zz 83 REAL(rstd) :: p(iim*jjm,llm+1) 84 REAL(rstd) :: Tave,Tp 85 REAL(rstd) :: T(iim*jjm,llm) 86 REAL(rstd) :: vt 87 REAL(rstd) :: d1,d2,d 88 REAL(rstd) :: ulon(3*iim*jjm,llm) 89 REAL(rstd) :: ulat(3*iim*jjm,llm) 90 REAL(rstd) :: lon,lat 91 REAL(rstd) :: pks(iim*jjm) 92 REAL(rstd) :: pk(iim*jjm,llm) 28 SUBROUTINE getin_etat0 29 dcmip5_testcase=1 30 CALL getin("dcmip5_testcase",dcmip5_testcase) 31 END SUBROUTINE getin_etat0 93 32 33 SUBROUTINE compute_etat0(ngrid,lon,lat, phis, ps, Temp, ulon, ulat, q) 34 USE disvert_mod 35 IMPLICIT NONE 36 INTEGER, INTENT(IN) :: ngrid 37 REAL(rstd),INTENT(IN) :: lon(ngrid) 38 REAL(rstd),INTENT(IN) :: lat(ngrid) 39 REAL(rstd),INTENT(OUT) :: ps(ngrid) 40 REAL(rstd),INTENT(OUT) :: phis(ngrid) 41 REAL(rstd),INTENT(OUT) :: Temp(ngrid,llm) 42 REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm) 43 REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm) 44 REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot) 94 45 95 latc=Pi/1896 lonc=Pi46 INTEGER :: l,ij 47 REAL(rstd) :: r, d, d1,d2, zz, Tave, Tp, vt, aa, bb 97 48 98 Tv0=T0*(1+0.608*q0)99 Tvt=Tv0-Gamma*zt100 101 testcase=1102 CALL getin("dcmip5_testcase",testcase)103 104 49 phis(:)=0. 105 106 DO j=jj_begin,jj_end 107 DO i=ii_begin,ii_end 108 ij=(j-1)*iim+i 109 CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 110 r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) 111 ps(ij)=pb-DeltaP*exp(-(r/rp)**1.5) 112 ENDDO 113 ENDDO 114 115 CALL compute_pression(ps,p,0) 50 51 DO ij=1,ngrid 52 r=radius*arc(lonc,latc, lon(ij),lat(ij)) 53 ps(ij)=pb-DeltaP*exp(-(r/rp)**1.5) 54 END DO 116 55 117 56 DO l=1,llm 118 DO j=jj_begin,jj_end 119 DO i=ii_begin,ii_end 120 ij=(j-1)*iim+i 121 CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 122 r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) 123 CALL compute_z(0.5*(p(ij,l)+p(ij,l+1)),ps(ij),r,z(ij,l)) 124 ENDDO 125 ENDDO 126 ENDDO 127 out_i=z 128 129 DO l=1,llm 130 DO j=jj_begin-1,jj_end+1 131 DO i=ii_begin-1,ii_end+1 132 ij=(j-1)*iim+i 133 IF (z(ij,l)<=zt) THEN 134 q(ij,l,1)=q0*exp(-z(ij,l)/zq1)*exp(-(z(ij,l)/zq2)**2) 57 aa=.5*(ap(l)+ap(l+1)) 58 bb=.5*(bp(l)+bp(l+1)) 59 DO ij=1,ngrid 60 r=radius*arc(lonc,latc, lon(ij),lat(ij)) 61 CALL compute_z(aa+bb*ps(ij),ps(ij),r,zz) 62 IF (zz<=zt) THEN 63 q(ij,l,1)=q0*exp(-zz/zq1)*exp(-(zz/zq2)**2) 64 Tave=Tv0-Gamma*zz 65 Tp=(Tv0-Gamma*zz) * ( ( 1 + 2*Rd*(Tv0-Gamma*zz)*zz / & 66 (g*zp**2*(1-pb/Deltap*exp((r/rp)**1.5+(zz/zp)**2))))**(-1) - 1) 67 vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 - 1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd & 68 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2) & 69 - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 ))) 135 70 ELSE 136 q(ij,l,1)=qt 71 q(ij,l,1)=qt 72 Tave=Tvt 73 Tp=0 74 vt=0 137 75 ENDIF 138 ENDDO 139 ENDDO 140 ENDDO 76 Temp(ij,l)=Tave+Tp 141 77 142 143 DO l=1,llm 144 DO j=jj_begin-1,jj_end+1 145 DO i=ii_begin-1,ii_end+1 146 ij=(j-1)*iim+i 147 CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 148 r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) 149 IF (z(ij,l)<=zt) THEN 150 Tave=Tv0-Gamma*z(ij,l) 151 Tp=(Tv0-Gamma*z(ij,l)) * ( ( 1 + 2*Rd*(Tv0-Gamma*z(ij,l))*z(ij,l) / & 152 (g*zp**2*(1-pb/Deltap*exp((r/rp)**1.5+(z(ij,l)/zp)**2))))**(-1) - 1) 153 ELSE 154 Tave=Tvt 155 Tp=0 156 ENDIF 157 T(ij,l)=Tave+Tp 158 out_i(ij,l)=Tave+Tp 159 ENDDO 160 ENDDO 161 ENDDO 162 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 168 CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat) 169 zz=0.5*(z(ij,l)+z(ij+t_right,l)) 170 r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) 171 IF (zz<=zt) THEN 172 vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 - 1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd & 173 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2) & 174 - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 ))) 175 ELSE 176 vt = 0 177 ENDIF 178 d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc) 179 d2=cos(latc)*sin(lon-lonc) 78 d1=sin(latc)*cos(lat(ij))-cos(latc)*sin(lat(ij))*cos(lon(ij)-lonc) 79 d2=cos(latc)*sin(lon(ij)-lonc) 180 80 d=MAX(1e-25,sqrt(d1**2+d2**2)) 181 81 ulon(ij+u_right,l)=vt*d1/d 182 82 ulat(ij+u_right,l)=vt*d2/d 183 184 185 186 CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat) 187 zz=0.5*(z(ij,l)+z(ij+t_lup,l)) 188 r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) 189 IF (zz<=zt) THEN 190 vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 - 1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd & 191 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2) & 192 - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 ))) 193 194 ELSE 195 vt = 0 196 ENDIF 197 d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc) 198 d2=cos(latc)*sin(lon-lonc) 199 d=MAX(1e-25,sqrt(d1**2+d2**2)) 200 ulon(ij+u_lup,l)=vt*d1/d 201 ulat(ij+u_lup,l)=vt*d2/d 202 203 204 205 206 CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat) 207 zz=0.5*(z(ij,l)+z(ij+t_ldown,l)) 208 r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) 209 IF (zz<=zt) THEN 210 vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 - 1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd & 211 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2) & 212 - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 ))) 213 214 ELSE 215 vt = 0 216 ENDIF 217 d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc) 218 d2=cos(latc)*sin(lon-lonc) 219 d=MAX(1e-25,sqrt(d1**2+d2**2)) 220 ulon(ij+u_ldown,l)=vt*d1/d 221 ulat(ij+u_ldown,l)=vt*d2/d 222 223 224 225 CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 226 zz=z(ij,l) 227 r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) 228 IF (zz<=zt) THEN 229 vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 - 1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd & 230 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2) & 231 - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 ))) 232 ELSE 233 vt = 0 234 ENDIF 235 d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc) 236 d2=cos(latc)*sin(lon-lonc) 237 d=MAX(1e-25,sqrt(d1**2+d2**2)) 238 ! ulon(ij+u_ldown,l)=vt*d1/d 239 ! ulat(ij+u_ldown,l)=vt*d2/d 240 ! out_i(ij,l)=sqrt((vt*d1/d)**2+(vt*d2/d)**2) 241 ENDDO 242 ENDDO 83 END DO 243 84 ENDDO 244 245 CALL compute_wind_perp_from_lonlat_compound(ulon, ulat, ue) 246 CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0) 247 248 END SUBROUTINE compute_etat0_dcmip5 85 END SUBROUTINE compute_etat0 249 86 250 87 SUBROUTINE compute_z(pmodel,ps,r,z) … … 257 94 REAL(rstd),INTENT(OUT) :: z 258 95 259 REAL(rstd) :: Tv0,Tvt,pt 260 REAL(rstd) :: pave, pp, dpdz 96 REAL(rstd) :: pt, pave, pp, dpdz 261 97 REAL(rstd) :: znp1 262 98 REAL(rstd) :: epsilon 263 99 REAL(rstd) :: p 264 INTEGER :: n 265 266 Tv0=T0/(1+0.608*q0) 267 Tvt=Tv0-Gamma*zt 100 INTEGER :: n 101 268 102 pt=pb*(Tvt/Tv0)**(g/(Rd*Gamma)) 269 103 … … 274 108 ENDIF 275 109 276 ! PRINT*,pmodel,pt,z,zt277 110 IF (z>zt .OR. r>1000000) return 278 111 … … 297 130 298 131 epsilon=ABS( (znp1-z)/znp1 ) 299 ! PRINT *,"epsilon",epsilon,r,z,znp1300 ! PRINT *,"----",n,pmodel,ps,pave,dpdz301 132 z=znp1 302 133 n=n+1 -
codes/icosagcm/trunk/src/etat0_jablonowsky06.f90
r201 r203 11 11 REAL(rstd),PARAMETER :: Gamma=0.005 12 12 REAL(rstd),PARAMETER :: up0=1 13 PUBLIC test_etat0_jablonowsky06, etat0, compute_etat0_jablonowsky06, compute_etat0_new 13 14 PUBLIC compute_etat0 15 14 16 CONTAINS 15 17 16 SUBROUTINE test_etat0_jablonowsky06 17 USE icosa 18 USE kinetic_mod 19 USE pression_mod 20 USE exner_mod 21 USE geopotential_mod 22 USE vorticity_mod 23 IMPLICIT NONE 24 TYPE(t_field),POINTER,SAVE :: f_ps(:) 25 TYPE(t_field),POINTER,SAVE :: f_phis(:) 26 TYPE(t_field),POINTER,SAVE :: f_theta_rhodz(:) 27 TYPE(t_field),POINTER,SAVE :: f_u(:) 28 TYPE(t_field),POINTER,SAVE :: f_Ki(:) 29 TYPE(t_field),POINTER,SAVE :: f_temp(:) 30 TYPE(t_field),POINTER,SAVE :: f_p(:) 31 TYPE(t_field),POINTER,SAVE :: f_pks(:) 32 TYPE(t_field),POINTER,SAVE :: f_pk(:) 33 TYPE(t_field),POINTER,SAVE :: f_phi(:) 34 TYPE(t_field),POINTER,SAVE :: f_vort(:) 35 TYPE(t_field),POINTER,SAVE :: f_q(:) 36 37 REAL(rstd),POINTER :: Ki(:,:) 38 REAL(rstd),POINTER :: temp(:) 39 INTEGER :: ind 40 41 42 CALL allocate_field(f_ps,field_t,type_real) 43 CALL allocate_field(f_phis,field_t,type_real) 44 CALL allocate_field(f_theta_rhodz,field_t,type_real,llm) 45 CALL allocate_field(f_u,field_u,type_real,llm) 46 CALL allocate_field(f_p,field_t,type_real,llm+1) 47 CALL allocate_field(f_Ki,field_t,type_real,llm) 48 CALL allocate_field(f_pks,field_t,type_real) 49 CALL allocate_field(f_pk,field_t,type_real,llm) 50 CALL allocate_field(f_phi,field_t,type_real,llm) 51 CALL allocate_field(f_temp,field_t,type_real) 52 CALL allocate_field(f_vort,field_z,type_real,llm) 53 54 CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 55 56 CALL kinetic(f_u,f_Ki) 57 CALL vorticity(f_u,f_vort) 58 59 CALL pression(f_ps,f_p) 60 CALL exner(f_ps,f_p,f_pks,f_pk) 61 CALL geopotential(f_phis,f_pks,f_pk,f_theta_rhodz,f_phi) 62 63 CALL writefield('ps',f_ps) 64 CALL writefield('phis',f_phis) 65 CALL writefield('theta',f_theta_rhodz) 66 CALL writefield('f_phi',f_phi) 67 68 CALL writefield('Ki',f_Ki) 69 CALL writefield('vort',f_vort) 70 71 END SUBROUTINE test_etat0_jablonowsky06 72 73 74 SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 75 USE icosa 76 IMPLICIT NONE 77 TYPE(t_field),POINTER :: f_ps(:) 78 TYPE(t_field),POINTER :: f_phis(:) 79 TYPE(t_field),POINTER :: f_theta_rhodz(:) 80 TYPE(t_field),POINTER :: f_u(:) 81 TYPE(t_field),POINTER :: f_q(:) 82 83 REAL(rstd),POINTER :: ps(:) 84 REAL(rstd),POINTER :: phis(:) 85 REAL(rstd),POINTER :: theta_rhodz(:,:) 86 REAL(rstd),POINTER :: u(:,:) 87 REAL(rstd),POINTER :: q(:,:,:) 88 INTEGER :: ind 89 90 DO ind=1,ndomain 91 IF (.NOT. assigned_domain(ind)) CYCLE 92 CALL swap_dimensions(ind) 93 CALL swap_geometry(ind) 94 ps=f_ps(ind) 95 phis=f_phis(ind) 96 theta_rhodz=f_theta_rhodz(ind) 97 u=f_u(ind) 98 q=f_q(ind) 99 q=0 100 CALL compute_etat0_jablonowsky06(ps, phis, theta_rhodz, u) 101 ENDDO 102 103 END SUBROUTINE etat0 104 105 SUBROUTINE compute_etat0_jablonowsky06(ps, phis, theta_rhodz, u) 106 USE icosa 107 USE disvert_mod 108 USE pression_mod 109 USE exner_mod 110 USE geopotential_mod 111 USE theta2theta_rhodz_mod 112 IMPLICIT NONE 113 REAL(rstd),INTENT(OUT) :: ps(iim*jjm) 114 REAL(rstd),INTENT(OUT) :: phis(iim*jjm) 115 REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) 116 REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm) 117 118 INTEGER :: i,j,l,ij 119 REAL(rstd) :: theta(iim*jjm,llm) 120 REAL(rstd) :: eta(llm) 121 REAL(rstd) :: etav(llm) 122 REAL(rstd) :: etas, etavs 123 REAL(rstd) :: lon,lat 124 REAL(rstd) :: ulon(3) 125 REAL(rstd) :: ep(3), norm_ep 126 REAL(rstd) :: Tave, T 127 REAL(rstd) :: phis_ave 128 REAL(rstd) :: V0(3) 129 REAL(rstd) :: r2 130 REAL(rstd) :: utot 131 REAL(rstd) :: lonx,latx 132 133 DO l=1,llm 134 eta(l)= 0.5 *( ap(l)/ps0+bp(l) + ap(l+1)/ps0+bp(l+1) ) 135 etav(l)=(eta(l)-eta0)*Pi/2 136 ENDDO 137 etas=ap(1)+bp(1) 138 etavs=(etas-eta0)*Pi/2 139 140 DO j=jj_begin,jj_end 141 DO i=ii_begin,ii_end 142 ij=(j-1)*iim+i 143 ps(ij)=ps0 144 ENDDO 145 ENDDO 146 147 148 CALL lonlat2xyz(Pi/9,2*Pi/9,V0) 149 150 u(:,:)=1e10 151 DO l=1,llm 152 DO j=jj_begin,jj_end 153 DO i=ii_begin,ii_end 154 ij=(j-1)*iim+i 155 156 CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon,lat) 157 CALL cross_product2(V0,xyz_e(ij+u_right,:)/radius,ep) 158 r2=(asin(sqrt(sum(ep*ep))))**2 159 utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01) 160 161 ulon(1) = -sin(lon) * utot 162 ulon(2) = cos(lon) * utot 163 ulon(3) = 0 * utot 164 165 166 CALL cross_product2(xyz_v(ij+z_rdown,:)/radius,xyz_v(ij+z_rup,:)/radius,ep) 167 norm_ep=sqrt(sum(ep(:)**2)) 168 IF (norm_ep>1e-30) THEN 169 ep=-ep/norm_ep*ne(ij,right) 170 u(ij+u_right,l)=sum(ep(:)*ulon(:)) 171 ENDIF 172 173 174 175 CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon,lat) 176 CALL cross_product2(V0,xyz_e(ij+u_lup,:)/radius,ep) 177 r2=(asin(sqrt(sum(ep*ep))))**2 178 utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01) 179 ulon(1) = -sin(lon) * utot 180 ulon(2) = cos(lon) * utot 181 ulon(3) = 0 * utot 182 183 184 CALL cross_product2(xyz_v(ij+z_up,:)/radius,xyz_v(ij+z_lup,:)/radius,ep) 185 norm_ep=sqrt(sum(ep(:)**2)) 186 IF (norm_ep>1e-30) THEN 187 ep=-ep/norm_ep*ne(ij,lup) 188 u(ij+u_lup,l)=sum(ep(:)*ulon(:)) 189 ENDIF 190 191 192 193 194 195 196 CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon,lat) 197 CALL cross_product2(V0,xyz_e(ij+u_ldown,:)/radius,ep) 198 r2=(asin(sqrt(sum(ep*ep))))**2 199 utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01) 200 ulon(1) = -sin(lon) * utot 201 ulon(2) = cos(lon) * utot 202 ulon(3) = 0 * utot 203 204 205 CALL cross_product2(xyz_v(ij+z_ldown,:)/radius,xyz_v(ij+z_down,:)/radius,ep) 206 norm_ep=sqrt(sum(ep(:)**2)) 207 IF (norm_ep>1e-30) THEN 208 ep=-ep/norm_ep*ne(ij,ldown) 209 u(ij+u_ldown,l)=sum(ep(:)*ulon(:)) 210 ENDIF 211 212 ENDDO 213 ENDDO 214 ENDDO 215 216 217 DO l=1,llm 218 Tave=T0*eta(l)**(Rd*Gamma/g) 219 IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5 220 DO j=jj_begin,jj_end 221 DO i=ii_begin,ii_end 222 ij=(j-1)*iim+i 223 CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) 224 225 T=Tave+ 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5 & 226 * ( (-2*sin(lat)**6*(cos(lat)**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5 & 227 + (8./5*cos(lat)**3*(sin(lat)**2+2./3)-Pi/4)*radius*Omega) 228 229 theta(ij,l)=T*eta(l)**(-kappa) 230 231 ENDDO 232 ENDDO 233 ENDDO 234 235 236 phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g)) 237 DO j=jj_begin,jj_end 238 DO i=ii_begin,ii_end 239 ij=(j-1)*iim+i 240 CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) 241 phis(ij)=phis_ave+u0*cos(etavs)**1.5*( (-2*sin(lat)**6 * (cos(lat)**2+1./3) + 10./63 )*u0*cos(etavs)**1.5 & 242 +(8./5*cos(lat)**3 * (sin(lat)**2 + 2./3) - Pi/4)*radius*Omega ) 243 ENDDO 244 ENDDO 245 246 CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0) 247 248 END SUBROUTINE compute_etat0_jablonowsky06 249 250 SUBROUTINE compute_etat0_new(ngrid,lat,lon, phis, ps, temp, ulon, ulat) 18 SUBROUTINE compute_etat0(ngrid,lon,lat, phis, ps, temp, ulon, ulat) 251 19 USE disvert_mod 252 20 IMPLICIT NONE … … 283 51 IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5 284 52 DO ij=1,ngrid 285 r2 = arc(Pi/9 ,2*Pi/9, lon(ij),lat(ij))**253 r2 = arc(Pi/9.,2.*Pi/9., lon(ij),lat(ij))**2 286 54 temp(ij,l) = Tave + 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5 & 287 55 * ( (-2*sin(lat(ij))**6*(cos(lat(ij))**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5 & … … 293 61 294 62 295 END SUBROUTINE compute_etat0 _new63 END SUBROUTINE compute_etat0 296 64 297 65 END MODULE etat0_jablonowsky06_mod -
codes/icosagcm/trunk/src/vector.f90
r201 r203 48 48 END SUBROUTINE cross_product2 49 49 50 FUNCTION arc(lon1,lat1, lon2,lat2) 51 REAL(rstd) :: lon1,lat1, lon2,lat2, arc2, & 52 x1,x2,x3,y1,y2,y3,z1,z2,z3 53 x1=cos(lat1)*cos(lon1) 54 x2=cos(lat2)*cos(lon2) 55 y1=cos(lat1)*sin(lon1) 56 y2=cos(lat2)*sin(lon2) 57 z1=sin(lat1) 58 z2=sin(lat2) 59 x3=y1*z2-y2*z1 60 y3=z1*x2-z2*x1 61 z3=x1*y2-x2*y1 62 arc=ASIN(SQRT(x3*x3+y3*y3+z3*z3)) 50 FUNCTION arc(lon,lat, lonc,latc) 51 REAL(rstd) :: lon,lat, lonc,latc, arc 52 arc=ACOS(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) 63 53 END FUNCTION arc 64 54
Note: See TracChangeset
for help on using the changeset viewer.
