Changeset 153 for codes/icosagcm/trunk
- Timestamp:
- 05/17/13 18:21:50 (11 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/geometry.f90
r151 r153 55 55 INTEGER, PARAMETER :: ne_ldown=1 56 56 INTEGER, PARAMETER :: ne_rdown=-1 57 57 58 REAL(rstd),private :: schmidt_factor, schmidt_lon, schmidt_lat 59 58 60 CONTAINS 59 61 … … 115 117 !$OMP BARRIER 116 118 END SUBROUTINE swap_geometry 117 119 120 SUBROUTINE update_circumcenters 121 USE domain_mod 122 USE dimensions 123 USE spherical_geom_mod 124 USE vector 125 USE transfert_mod 126 127 IMPLICIT NONE 128 REAL(rstd) :: x1(3),x2(3) 129 REAL(rstd) :: vect(3,6) 130 REAL(rstd) :: centr(3) 131 132 INTEGER :: ind,i,j,n,k 133 134 CALL transfert_request(geom%xyz_i,req_i1) 135 136 DO ind=1,ndomain 137 CALL swap_dimensions(ind) 138 CALL swap_geometry(ind) 139 DO j=jj_begin,jj_end 140 DO i=ii_begin,ii_end 141 n=(j-1)*iim+i 142 DO k=0,5 143 x1(:) = xyz_i(n+t_pos(k+1),:) 144 x2(:) = xyz_i(n+t_pos(MOD(k+1,6)+1),:) 145 if (norm(x1-x2)<1e-16) x2(:) = xyz_i(n+t_pos(MOD(k+2,6)+1),:) 146 CALL circumcenter(xyz_i(n,:), x1, x2, xyz_v(n+z_pos(k+1),:)) 147 ENDDO 148 ENDDO 149 ENDDO 150 ENDDO 151 152 END SUBROUTINE update_circumcenters 153 118 154 SUBROUTINE optimize_geometry 119 155 USE metric … … 149 185 ENDDO 150 186 151 CALL transfert_request(geom%xyz_i,req_i1) 152 153 DO ind=1,ndomain 154 CALL swap_dimensions(ind) 155 CALL swap_geometry(ind) 156 DO j=jj_begin,jj_end 157 DO i=ii_begin,ii_end 158 n=(j-1)*iim+i 159 DO k=0,5 160 x1(:) = xyz_i(n+t_pos(k+1),:) 161 x2(:) = xyz_i(n+t_pos(MOD(k+1,6)+1),:) 162 if (norm(x1-x2)<1e-16) x2(:) = xyz_i(n+t_pos(MOD(k+2,6)+1),:) 163 CALL circumcenter(xyz_i(n,:), x1, x2, xyz_v(n+z_pos(k+1),:)) 164 ENDDO 165 ENDDO 166 ENDDO 167 ENDDO 168 187 CALL update_circumcenters 188 169 189 DO ind=1,ndomain 170 190 d=>domain(ind) … … 221 241 PRINT *,"it = ",it," diff centroid circumcenter ",sum 222 242 ENDIF 243 244 CALL update_circumcenters 245 246 ENDDO 247 248 END SUBROUTINE optimize_geometry 249 250 SUBROUTINE schmidt_transform(xyz,cc) 251 ! Based on formula (12) from Guo & Drake, JCP 2005 252 USE spherical_geom_mod 253 IMPLICIT NONE 254 REAL(rstd),INTENT(INOUT) :: xyz(3) 255 REAL(rstd), INTENT(IN) :: cc ! stretching factor>0 256 REAL(rstd) :: lat,lon,mu 257 CALL xyz2lonlat_relative(xyz,schmidt_lon,schmidt_lat, lon,lat) 258 mu = sin(lat) 259 mu = ((cc-1)+mu*(cc+1)) / ((cc+1)+mu*(cc-1)) 260 lat = asin(mu) 261 CALL lonlat2xyz_relative(lon,lat,schmidt_lon, schmidt_lat, xyz) 262 END SUBROUTINE schmidt_transform 223 263 224 CALL transfert_request(geom%xyz_i,req_i1)225 226 DO ind=1,ndomain227 CALL swap_dimensions(ind)228 CALL swap_geometry(ind)229 DO j=jj_begin,jj_end230 DO i=ii_begin,ii_end231 n=(j-1)*iim+i232 DO k=0,5233 x1(:) = xyz_i(n+t_pos(k+1),:)234 x2(:) = xyz_i(n+t_pos(MOD(k+1,6)+1),:)235 if (norm(x1-x2)<1e-16) x2(:) = xyz_i(n+t_pos(MOD(k+2,6)+1),:)236 CALL circumcenter(xyz_i(n,:), x1, x2, xyz_v(n+z_pos(k+1),:))237 ENDDO238 ENDDO239 ENDDO240 ENDDO241 242 243 ENDDO244 245 END SUBROUTINE optimize_geometry246 247 248 264 SUBROUTINE set_geometry 249 265 USE metric … … 253 269 USE dimensions 254 270 USE transfert_mod 271 USE ioipsl 255 272 IMPLICIT NONE 256 273 … … 259 276 REAL(rstd) :: vect(3,6) 260 277 REAL(rstd) :: centr(3) 261 REAL(rstd) :: vet(3),vep(3) 278 REAL(rstd) :: vet(3),vep(3), vertex(3) 262 279 INTEGER :: ind,i,j,k,n 263 280 TYPE(t_domain),POINTER :: d … … 270 287 271 288 CALL optimize_geometry 272 289 290 ! Schmidt transform 291 schmidt_factor = 1. 292 CALL getin('schmidt_factor', schmidt_factor) 293 294 schmidt_lon = 0. 295 CALL getin('schmidt_lon', schmidt_lon) 296 schmidt_lon = schmidt_lon * pi/180. 297 298 schmidt_lat = 45. 299 CALL getin('schmidt_lat', schmidt_lat) 300 schmidt_lat = schmidt_lat * pi/180. 301 302 DO ind=1,ndomain 303 CALL swap_dimensions(ind) 304 CALL swap_geometry(ind) 305 DO j=jj_begin-1,jj_end+1 306 DO i=ii_begin-1,ii_end+1 307 n=(j-1)*iim+i 308 vertex(:)=xyz_i(n,:) 309 CALL schmidt_transform(vertex, schmidt_factor**2.) 310 xyz_i(n,:)=vertex(:) 311 END DO 312 END DO 313 END DO 314 315 ! at this point the mesh is defined on the unit sphere by xyz_i 316 ! all other mesh quantities are now deduced from xyz_i 317 318 CALL update_circumcenters 319 320 ! DO ind=1,ndomain 321 ! d=>domain(ind) 322 ! CALL swap_dimensions(ind) 323 ! CALL swap_geometry(ind) 324 ! DO j=1,jjm 325 ! DO i=1,iim 326 ! n=(j-1)*iim+i 327 ! d%xyz(:,i,j)=xyz_i(n,:) 328 ! ENDDO 329 ! ENDDO 330 ! ENDDO 331 ! 332 ! CALL compute_boundary 333 273 334 DO ind=1,ndomain 274 335 d=>domain(ind) -
codes/icosagcm/trunk/src/spherical_geom.f90
r15 r153 34 34 35 35 lat=asin(xyzn(3)) 36 coslat=cos(lat) 37 IF (abs(coslat)<1e-15) THEN 38 lon=0. 39 ELSE 40 lon=acos(MAX(MIN((xyzn(1)/coslat),1.),-1.)) 41 IF (xyzn(2)/coslat<0) lon=-lon 42 ENDIF 43 36 lon=atan2(xyzn(2),xyzn(1)) 44 37 END SUBROUTINE xyz2lonlat 38 39 ! lat/lon with respect to a displaced pole (rotated basis) : 40 ! ( cos(lon0)*sin(lat0), sin(lon0)*sin(lat0), -cos(lat0)) 41 ! (-sin(lon0), cos(lon0), 0) 42 ! ( cos(lon0)*cos(lat0), sin(lon0)*cos(lat0), sin(lat0)) 43 44 SUBROUTINE lonlat2xyz_relative(lon,lat,lon0,lat0, xyz) 45 IMPLICIT NONE 46 REAL(rstd),INTENT(IN) :: lon0, lat0, lon,lat 47 REAL(rstd),INTENT(OUT) :: xyz(3) 48 REAL(rstd) :: xx,yy,zz 49 xx = cos(lon)*cos(lat) 50 yy = sin(lon)*cos(lat) 51 zz = sin(lat) 52 xyz(1) = cos(lon0)*(sin(lat0)*xx+cos(lat0)*zz)-sin(lon0)*yy 53 xyz(2) = sin(lon0)*(sin(lat0)*yy+cos(lat0)*zz)+cos(lon0)*yy 54 xyz(3) = sin(lat0)*zz-cos(lat0)*xx 55 END SUBROUTINE lonlat2xyz_relative 56 57 SUBROUTINE xyz2lonlat_relative(xyz,lon0,lat0, lon,lat) 58 IMPLICIT NONE 59 REAL(rstd),INTENT(IN) :: xyz(3), lon0, lat0 60 REAL(rstd),INTENT(OUT) :: lon,lat 61 REAL(rstd) :: xx,yy,zz 62 xx = sin(lat0)*(xyz(1)*cos(lon0)+xyz(2)*sin(lon0))-cos(lat0)*xyz(3) 63 yy = xyz(2)*cos(lon0)-xyz(1)*sin(lon0) 64 zz = cos(lat0)*(xyz(1)*cos(lon0)+xyz(2)*sin(lon0))+sin(lat0)*xyz(3) 65 lon = atan2(yy,xx) 66 lat = asin(zz) 67 END SUBROUTINE xyz2lonlat_relative 68 45 69 46 70
Note: See TracChangeset
for help on using the changeset viewer.