Changeset 15 for codes/icosagcm/trunk
- Timestamp:
- 07/09/12 15:23:38 (12 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 3 added
- 1 deleted
- 16 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/caldyn_gcm.f90
r12 r15 11 11 REAL(rstd),POINTER :: out_z(:,:) 12 12 13 INTEGER :: itau_out 14 13 15 CONTAINS 14 16 17 SUBROUTINE init_caldyn(dt) 18 USE IOIPSL 19 IMPLICIT NONE 20 REAL(rstd),INTENT(IN) :: dt 21 REAL(rstd) :: write_period 22 CALL allocate_caldyn 23 24 CALL getin('write_period',write_period) 25 26 itau_out=INT(write_period/dt) 27 28 END SUBROUTINE init_caldyn 29 15 30 SUBROUTINE allocate_caldyn 16 31 USE domain_mod … … 78 93 79 94 80 SUBROUTINE caldyn( f_phis, f_ps, f_theta_rhodz, f_u, f_dps, f_dtheta_rhodz, f_du)95 SUBROUTINE caldyn(it,f_phis, f_ps, f_theta_rhodz, f_u, f_dps, f_dtheta_rhodz, f_du) 81 96 USE domain_mod 82 97 USE dimensions … … 87 102 USE vorticity_mod 88 103 USE kinetic_mod 104 USE theta2theta_rhodz_mod 89 105 IMPLICIT NONE 106 INTEGER,INTENT(IN) :: it 90 107 TYPE(t_field),POINTER :: f_phis(:) 91 108 TYPE(t_field),POINTER :: f_ps(:) … … 104 121 REAL(rstd),POINTER :: du(:,:) 105 122 INTEGER :: ind 106 INTEGER,SAVE :: it=0 123 107 124 108 125 CALL transfert_request(f_phis,req_i1) … … 139 156 ! CALL kinetic(f_du,f_out) 140 157 141 IF (mod(it, 72)==0 ) THEN158 IF (mod(it,itau_out)==0 ) THEN 142 159 CALL writefield("ps",f_ps) 143 CALL writefield("dps",f_dps)160 ! CALL writefield("dps",f_dps) 144 161 ! CALL writefield("theta_rhodz",f_theta_rhodz) 145 CALL writefield("dtheta_rhodz",f_dtheta_rhodz) 146 ! CALL writefield("vort",f_out_z) 147 CALL writefield("theta",f_out) 162 ! CALL writefield("dtheta_rhodz",f_dtheta_rhodz) 163 CALL vorticity(f_u,f_out_z) 164 CALL writefield("vort",f_out_z) 165 ! CALL writefield("theta",f_out) 166 CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_out) ; 167 CALL writefield("T",f_out) 148 168 149 169 ! CALL writefield("out",f_out) … … 155 175 ENDIF 156 176 ! CALL check_mass_conservation(f_ps,f_dps) 157 it=it+1 177 158 178 END SUBROUTINE caldyn 159 179 -
codes/icosagcm/trunk/src/caldyn_sw.f90
r12 r15 170 170 CALL writefield("dh",f_dh) 171 171 CALL Compute_enstrophy 172 ! CALL Compute_PV 172 173 ENDIF 173 174 it=it+1 -
codes/icosagcm/trunk/src/domain.f90
r12 r15 190 190 ng2=d%neighbour(:,MOD(k+1,6),i,j) 191 191 IF (sqrt(sum((ng1-ng2)**2))<1e-16) ng2=d%neighbour(:,MOD(k+2,6),i,j) 192 CALL circ onscrit(d%xyz(:,i,j),ng1,ng2,d%vertex(:,k,i,j))192 CALL circumcenter(d%xyz(:,i,j),ng1,ng2,d%vertex(:,k,i,j)) 193 193 ENDDO 194 194 ENDDO -
codes/icosagcm/trunk/src/etat0_jablonowsky06.f90
r12 r15 140 140 REAL(rstd) :: r2 141 141 REAL(rstd) :: utot 142 REAL(rstd) :: lonx,latx 142 143 143 144 DO l=1,llm … … 157 158 158 159 CALL lonlat2xyz(Pi/9,2*Pi/9,V0) 160 159 161 u(:,:)=1e10 160 162 DO l=1,llm … … 165 167 CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon,lat) 166 168 CALL cross_product2(V0,xyz_e(ij+u_right,:)/radius,ep) 167 r2=(asin(s um(ep*ep)))**2169 r2=(asin(sqrt(sum(ep*ep))))**2 168 170 utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01) 169 171 … … 184 186 CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon,lat) 185 187 CALL cross_product2(V0,xyz_e(ij+u_lup,:)/radius,ep) 186 r2=(asin(s um(ep*ep)))**2188 r2=(asin(sqrt(sum(ep*ep))))**2 187 189 utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01) 188 190 ulon(1) = -sin(lon) * utot … … 205 207 CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon,lat) 206 208 CALL cross_product2(V0,xyz_e(ij+u_ldown,:)/radius,ep) 207 r2=(asin(s um(ep*ep)))**2209 r2=(asin(sqrt(sum(ep*ep))))**2 208 210 utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01) 209 211 ulon(1) = -sin(lon) * utot -
codes/icosagcm/trunk/src/geometry.f90
r12 r15 4 4 TYPE t_geometry 5 5 TYPE(t_field),POINTER :: xyz_i(:) 6 TYPE(t_field),POINTER :: centroid(:) 6 7 TYPE(t_field),POINTER :: xyz_e(:) 7 8 TYPE(t_field),POINTER :: xyz_v(:) 9 TYPE(t_field),POINTER :: ep_e(:) 10 TYPE(t_field),POINTER :: et_e(:) 11 TYPE(t_field),POINTER :: elon_i(:) 12 TYPE(t_field),POINTER :: elat_i(:) 13 TYPE(t_field),POINTER :: elon_e(:) 14 TYPE(t_field),POINTER :: elat_e(:) 8 15 TYPE(t_field),POINTER :: Ai(:) 9 16 TYPE(t_field),POINTER :: Av(:) … … 20 27 TYPE(t_geometry),TARGET :: geom 21 28 22 REAL(rstd),POINTER :: Ai(:) 23 REAL(rstd),POINTER :: xyz_i(:,:) 24 REAL(rstd),POINTER :: xyz_e(:,:) 25 REAL(rstd),POINTER :: xyz_v(:,:) 26 REAL(rstd),POINTER :: Av(:) 27 REAL(rstd),POINTER :: de(:) 28 REAL(rstd),POINTER :: le(:) 29 REAL(rstd),POINTER :: Riv(:,:) 30 REAL(rstd),POINTER :: Riv2(:,:) 31 INTEGER,POINTER :: ne(:,:) 32 REAL(rstd),POINTER :: Wee(:,:,:) 33 REAL(rstd),POINTER :: bi(:) 34 REAL(rstd),POINTER :: fv(:) 29 REAL(rstd),POINTER :: Ai(:) ! area of a cell 30 REAL(rstd),POINTER :: xyz_i(:,:) ! coordinate of the center of the cell (voronoi) 31 REAL(rstd),POINTER :: centroid(:,:) ! coordinate of the centroid of the cell 32 REAL(rstd),POINTER :: xyz_e(:,:) ! coordinate of a wind point on the cell on a edge 33 REAL(rstd),POINTER :: ep_e(:,:) ! perpendicular unit vector of a edge (outsider) 34 REAL(rstd),POINTER :: et_e(:,:) ! tangeantial unit vector of a edge 35 REAL(rstd),POINTER :: elon_i(:,:) ! unit longitude vector on the center 36 REAL(rstd),POINTER :: elat_i(:,:) ! unit latitude vector on the center 37 REAL(rstd),POINTER :: elon_e(:,:) ! unit longitude vector on a wind point 38 REAL(rstd),POINTER :: elat_e(:,:) ! unit latitude vector on a wind point 39 REAL(rstd),POINTER :: xyz_v(:,:) ! coordinate of a vertex (center of the dual mesh) 40 REAL(rstd),POINTER :: Av(:) ! area of dual mesk cell 41 REAL(rstd),POINTER :: de(:) ! distance from a neighbour == lenght of an edge of the dual mesh 42 REAL(rstd),POINTER :: le(:) ! lenght of a edge 43 REAL(rstd),POINTER :: Riv(:,:) ! weight 44 REAL(rstd),POINTER :: Riv2(:,:) ! weight 45 INTEGER,POINTER :: ne(:,:) ! convention for the way on the normal wind on an edge 46 REAL(rstd),POINTER :: Wee(:,:,:) ! weight 47 REAL(rstd),POINTER :: bi(:) ! orographie 48 REAL(rstd),POINTER :: fv(:) ! coriolis (evaluted on a vertex) 35 49 36 50 … … 44 58 CALL allocate_field(geom%Ai,field_t,type_real) 45 59 CALL allocate_field(geom%xyz_i,field_t,type_real,3) 60 CALL allocate_field(geom%centroid,field_t,type_real,3) 46 61 CALL allocate_field(geom%xyz_e,field_u,type_real,3) 62 CALL allocate_field(geom%ep_e,field_u,type_real,3) 63 CALL allocate_field(geom%et_e,field_u,type_real,3) 64 CALL allocate_field(geom%elon_i,field_t,type_real,3) 65 CALL allocate_field(geom%elat_i,field_t,type_real,3) 66 CALL allocate_field(geom%elon_e,field_u,type_real,3) 67 CALL allocate_field(geom%elat_e,field_u,type_real,3) 47 68 CALL allocate_field(geom%xyz_v,field_z,type_real,3) 48 69 CALL allocate_field(geom%de,field_u,type_real) … … 66 87 Ai=geom%Ai(ind) 67 88 xyz_i=geom%xyz_i(ind) 89 centroid=geom%centroid(ind) 68 90 xyz_e=geom%xyz_e(ind) 91 ep_e=geom%ep_e(ind) 92 et_e=geom%et_e(ind) 93 elon_i=geom%elon_i(ind) 94 elat_i=geom%elat_i(ind) 95 elon_e=geom%elon_e(ind) 96 elat_e=geom%elat_e(ind) 69 97 xyz_v=geom%xyz_v(ind) 70 98 de=geom%de(ind) … … 79 107 80 108 END SUBROUTINE swap_geometry 81 82 SUBROUTINE set_geometry109 110 SUBROUTINE optimize_geometry 83 111 USE metric 84 112 USE spherical_geom_mod 85 113 USE domain_mod 86 114 USE dimensions 115 USE transfert_mod 116 USE vector 117 IMPLICIT NONE 118 INTEGER,PARAMETER :: nb_it=3000 119 TYPE(t_domain),POINTER :: d 120 INTEGER :: ind,it,i,j,n,k 121 REAL(rstd) :: x1(3),x2(3) 122 REAL(rstd) :: vect(3,6) 123 REAL(rstd) :: centr(3) 124 REAL(rstd) :: sum 125 LOGICAL :: check 126 127 DO ind=1,ndomain 128 d=>domain(ind) 129 CALL swap_dimensions(ind) 130 CALL swap_geometry(ind) 131 DO j=jj_begin,jj_end 132 DO i=ii_begin,ii_end 133 n=(j-1)*iim+i 134 xyz_i(n,:)=d%xyz(:,i,j) 135 ENDDO 136 ENDDO 137 ENDDO 138 139 CALL transfert_request(geom%xyz_i,req_i1) 140 141 DO ind=1,ndomain 142 CALL swap_dimensions(ind) 143 CALL swap_geometry(ind) 144 DO j=jj_begin,jj_end 145 DO i=ii_begin,ii_end 146 n=(j-1)*iim+i 147 DO k=0,5 148 x1(:) = xyz_i(n+t_pos(k+1),:) 149 x2(:) = xyz_i(n+t_pos(MOD(k+1,6)+1),:) 150 if (norm(x1-x2)<1e-16) x2(:) = xyz_i(n+t_pos(MOD(k+2,6)+1),:) 151 CALL circumcenter(xyz_i(n,:), x1, x2, xyz_v(n+z_pos(k+1),:)) 152 ENDDO 153 ENDDO 154 ENDDO 155 ENDDO 156 157 DO ind=1,ndomain 158 d=>domain(ind) 159 CALL swap_dimensions(ind) 160 CALL swap_geometry(ind) 161 DO j=jj_begin,jj_end 162 DO i=ii_begin,ii_end 163 n=(j-1)*iim+i 164 DO k=0,5 165 x1(:) = xyz_v(n+z_pos(k+1),:) 166 x2(:) = d%vertex(:,k,i,j) 167 IF (norm(x1-x2)>1e-10) THEN 168 PRINT*,"vertex diff ",ind,i,j,k 169 PRINT*,x1 170 PRINT*,x2 171 ENDIF 172 ENDDO 173 ENDDO 174 ENDDO 175 ENDDO 176 177 178 DO it=1,nb_it 179 IF (MOD(it,100)==0) THEN 180 check=.TRUE. 181 ELSE 182 check=.FALSE. 183 ENDIF 184 185 sum=0 186 DO ind=1,ndomain 187 CALL swap_dimensions(ind) 188 CALL swap_geometry(ind) 189 DO j=jj_begin,jj_end 190 DO i=ii_begin,ii_end 191 n=(j-1)*iim+i 192 vect(:,1)=xyz_v(n+z_rup,:) 193 vect(:,2)=xyz_v(n+z_up,:) 194 vect(:,3)=xyz_v(n+z_lup,:) 195 vect(:,4)=xyz_v(n+z_ldown,:) 196 vect(:,5)=xyz_v(n+z_down,:) 197 vect(:,6)=xyz_v(n+z_rdown,:) 198 CALL compute_centroid(vect,6,centr) 199 IF (check) THEN 200 sum=MAX(sum,norm(xyz_i(n,:)-centr(:))) 201 ENDIF 202 xyz_i(n,:)=centr(:) 203 ENDDO 204 ENDDO 205 ! i=ii_begin ; j=jj_begin ; n=(j-1)*iim+i ; xyz_i(n,:)=domain(ind)%xyz(:,i,j) 206 ! i=ii_begin ; j=jj_end ; n=(j-1)*iim+i ; xyz_i(n,:)=domain(ind)%xyz(:,i,j) 207 ! i=ii_end ; j=jj_begin ; n=(j-1)*iim+i ; xyz_i(n,:)=domain(ind)%xyz(:,i,j) 208 ! PRINT *,"Pb ?? : " 209 ! PRINT *, xyz_i(n,:), domain(ind)%xyz(:,i,j), norm(xyz_i(n,:)- domain(ind)%xyz(:,i,j)) 210 ! i=ii_end ; j=jj_end ; n=(j-1)*iim+i ; xyz_i(n,:)=domain(ind)%xyz(:,i,j) 211 212 ENDDO 213 214 IF (check) THEN 215 ! sum=sum/(ndomain*ii_nb*jj_nb) 216 PRINT *,"it = ",it," diff centroid circumcenter ",sum 217 ENDIF 218 219 CALL transfert_request(geom%xyz_i,req_i1) 220 221 DO ind=1,ndomain 222 CALL swap_dimensions(ind) 223 CALL swap_geometry(ind) 224 DO j=jj_begin,jj_end 225 DO i=ii_begin,ii_end 226 n=(j-1)*iim+i 227 DO k=0,5 228 x1(:) = xyz_i(n+t_pos(k+1),:) 229 x2(:) = xyz_i(n+t_pos(MOD(k+1,6)+1),:) 230 if (norm(x1-x2)<1e-16) x2(:) = xyz_i(n+t_pos(MOD(k+2,6)+1),:) 231 CALL circumcenter(xyz_i(n,:), x1, x2, xyz_v(n+z_pos(k+1),:)) 232 ENDDO 233 ENDDO 234 ENDDO 235 ENDDO 236 237 238 ENDDO 239 240 END SUBROUTINE optimize_geometry 241 242 243 SUBROUTINE set_geometry 244 USE metric 245 USE vector 246 USE spherical_geom_mod 247 USE domain_mod 248 USE dimensions 249 USE transfert_mod 87 250 IMPLICIT NONE 88 251 89 252 REAL(rstd) :: surf(6) 90 253 REAL(rstd) :: surf_v(6) 254 REAL(rstd) :: vect(3,6) 255 REAL(rstd) :: centr(6) 256 REAL(rstd) :: vet(3),vep(3) 91 257 INTEGER :: ind,i,j,k,n 92 258 TYPE(t_domain),POINTER :: d 93 259 REAL(rstd) :: S 94 260 REAL(rstd) :: w(6) 95 REA l(rstd) :: lon,lat261 REAL(rstd) :: lon,lat 96 262 INTEGER :: ii_glo,jj_glo 97 263 REAL(rstd) :: S1,S2 264 98 265 266 CALL optimize_geometry 267 99 268 DO ind=1,ndomain 100 269 d=>domain(ind) … … 104 273 DO i=ii_begin-1,ii_end+1 105 274 n=(j-1)*iim+i 106 107 xyz_i(n,:)=d%xyz(:,i,j) 108 xyz_v(n+z_up,:)=d%vertex(:,vup-1,i,j) 109 xyz_v(n+z_down,:)=d%vertex(:,vdown-1,i,j) 110 275 276 DO k=0,5 277 ne(n,k+1)=d%ne(k,i,j) 278 ENDDO 279 280 ! xyz_i(n,:)=d%xyz(:,i,j) 281 ! xyz_v(n+z_up,:)=d%vertex(:,vup-1,i,j) 282 ! xyz_v(n+z_down,:)=d%vertex(:,vdown-1,i,j) 283 284 vect(:,1)=xyz_v(n+z_rup,:) 285 vect(:,2)=xyz_v(n+z_up,:) 286 vect(:,3)=xyz_v(n+z_lup,:) 287 vect(:,4)=xyz_v(n+z_ldown,:) 288 vect(:,5)=xyz_v(n+z_down,:) 289 vect(:,6)=xyz_v(n+z_rdown,:) 290 CALL compute_centroid(vect,6,centr) 291 centroid(n,:)=centr(:) 292 293 111 294 CALL xyz2lonlat(xyz_v(n+z_up,:),lon,lat) 112 295 fv(n+z_up)=2*sin(lat)*omega … … 116 299 bi(n)=0. 117 300 118 CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),de(n+u_right)) 119 CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),de(n+u_lup)) 120 CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),de(n+u_ldown)) 121 122 CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),0.5,xyz_e(n+u_right,:)) 123 CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),0.5,xyz_e(n+u_lup,:)) 124 CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),0.5,xyz_e(n+u_ldown,:)) 125 126 CALL dist_cart(d%vertex(:,vrdown-1,i,j),d%vertex(:,vrup-1,i,j),le(n+u_right)) 127 CALL dist_cart(d%vertex(:,vup-1,i,j),d%vertex(:,vlup-1,i,j),le(n+u_lup)) 128 CALL dist_cart(d%vertex(:,vldown-1,i,j),d%vertex(:,vdown-1,i,j),le(n+u_ldown)) 129 301 ! CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),de(n+u_right)) 302 ! CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),de(n+u_lup)) 303 ! CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),de(n+u_ldown)) 304 305 CALL dist_cart(xyz_i(n,:),xyz_i(n+t_right,:),de(n+u_right)) 306 CALL dist_cart(xyz_i(n,:),xyz_i(n+t_lup,:),de(n+u_lup)) 307 CALL dist_cart(xyz_i(n,:),xyz_i(n+t_ldown,:),de(n+u_ldown)) 308 309 ! CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),0.5,xyz_e(n+u_right,:)) 310 ! CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),0.5,xyz_e(n+u_lup,:)) 311 ! CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),0.5,xyz_e(n+u_ldown,:)) 312 313 CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_right,:),0.5,xyz_e(n+u_right,:)) 314 CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_lup,:),0.5,xyz_e(n+u_lup,:)) 315 CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_ldown,:),0.5,xyz_e(n+u_ldown,:)) 316 317 ! CALL dist_cart(d%vertex(:,vrdown-1,i,j),d%vertex(:,vrup-1,i,j),le(n+u_right)) 318 ! CALL dist_cart(d%vertex(:,vup-1,i,j),d%vertex(:,vlup-1,i,j),le(n+u_lup)) 319 ! CALL dist_cart(d%vertex(:,vldown-1,i,j),d%vertex(:,vdown-1,i,j),le(n+u_ldown)) 320 321 CALL dist_cart(xyz_v(n+z_rdown,:), xyz_v(n+z_rup,:),le(n+u_right)) 322 CALL dist_cart(xyz_v(n+z_up,:), xyz_v(n+z_lup,:),le(n+u_lup)) 323 CALL dist_cart(xyz_v(n+z_ldown,:), xyz_v(n+z_down,:),le(n+u_ldown)) 130 324 131 325 Ai(n)=0 132 326 DO k=0,5 133 CALL surf_triangle(d%xyz(:,i,j),d%neighbour(:,k,i,j),d%neighbour(:,MOD((k+1+6),6),i,j),surf_v(k+1)) 134 CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,MOD((k-1+6),6),i,j),d%vertex(:,k,i,j),surf(k+1)) 135 ne(n,k+1)=d%ne(k,i,j) 327 ! CALL surf_triangle(d%xyz(:,i,j),d%neighbour(:,k,i,j),d%neighbour(:,MOD((k+1+6),6),i,j),surf_v(k+1)) 328 ! CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,MOD((k-1+6),6),i,j),d%vertex(:,k,i,j),surf(k+1)) 329 CALL surf_triangle(xyz_i(n,:),xyz_i(n+t_pos(k+1),:),xyz_i(n+t_pos(MOD((k+1+6),6)+1),:),surf_v(k+1)) 330 CALL surf_triangle(xyz_i(n,:),xyz_v(n+z_pos(MOD((k-1+6),6)+1),:),xyz_v(n+z_pos(k+1),:),surf(k+1)) 136 331 Ai(n)=Ai(n)+surf(k+1) 137 ENDDO 138 139 332 IF (i==ii_end .AND. j==jj_begin) THEN 333 IF (Ai(n)<1e20) THEN 334 ELSE 335 PRINT *,"PB !!",Ai(n),k,surf(k+1) 336 PRINT*,xyz_i(n,:),xyz_v(n+z_pos(MOD((k-1+6),6)+1),:),xyz_v(n+z_pos(k+1),:) 337 ENDIF 338 ENDIF 339 ENDDO 340 341 342 CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),vep) 343 IF (norm(vep)>1e-30) THEN 344 vep(:)=vep(:)/norm(vep) 345 CALL cross_product2(vep,xyz_e(n+u_right,:),vet) 346 vet(:)=vet(:)/norm(vet) 347 ep_e(n+u_right,:)=vep(:)*ne(n,right) 348 et_e(n+u_right,:)=vet(:)*ne(n,right) 349 ENDIF 350 351 CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),vep) 352 IF (norm(vep)>1e-30) THEN 353 vep(:)=vep(:)/norm(vep) 354 CALL cross_product2(vep,xyz_e(n+u_lup,:),vet) 355 vet(:)=vet(:)/norm(vet) 356 ep_e(n+u_lup,:)=vep(:)*ne(n,lup) 357 et_e(n+u_lup,:)=vet(:)*ne(n,lup) 358 ENDIF 359 360 CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),vep) 361 IF (norm(vep)>1e-30) THEN 362 vep(:)=vep(:)/norm(vep) 363 CALL cross_product2(vep,xyz_e(n+u_ldown,:),vet) 364 vet(:)=vet(:)/norm(vet) 365 ep_e(n+u_ldown,:)=vep(:)*ne(n,ldown) 366 et_e(n+u_ldown,:)=vet(:)*ne(n,ldown) 367 ENDIF 368 369 CALL xyz2lonlat(xyz_i(n,:),lon,lat) 370 elon_i(n,1) = -sin(lon) 371 elon_i(n,2) = cos(lon) 372 elon_i(n,3) = 0 373 elat_i(n,1) = -cos(lon)*sin(lat) 374 elat_i(n,2) = -sin(lon)*sin(lat) 375 elat_i(n,3) = cos(lat) 376 377 378 CALL xyz2lonlat(xyz_e(n+u_right,:),lon,lat) 379 elon_e(n+u_right,1) = -sin(lon) 380 elon_e(n+u_right,2) = cos(lon) 381 elon_e(n+u_right,3) = 0 382 elat_e(n+u_right,1) = -cos(lon)*sin(lat) 383 elat_e(n+u_right,2) = -sin(lon)*sin(lat) 384 elat_e(n+u_right,3) = cos(lat) 385 386 CALL xyz2lonlat(xyz_e(n+u_lup,:),lon,lat) 387 elon_e(n+u_lup,1) = -sin(lon) 388 elon_e(n+u_lup,2) = cos(lon) 389 elon_e(n+u_lup,3) = 0 390 elat_e(n+u_lup,1) = -cos(lon)*sin(lat) 391 elat_e(n+u_lup,2) = -sin(lon)*sin(lat) 392 elat_e(n+u_lup,3) = cos(lat) 393 394 CALL xyz2lonlat(xyz_e(n+u_ldown,:),lon,lat) 395 elon_e(n+u_ldown,1) = -sin(lon) 396 elon_e(n+u_ldown,2) = cos(lon) 397 elon_e(n+u_ldown,3) = 0 398 elat_e(n+u_ldown,1) = -cos(lon)*sin(lat) 399 elat_e(n+u_ldown,2) = -sin(lon)*sin(lat) 400 elat_e(n+u_ldown,3) = cos(lat) 401 402 140 403 DO k=0,5 141 CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,k,i,j),S1) 142 CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,MOD(k+1+6,6),i,j),S2) 404 ! CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,k,i,j),S1) 405 ! CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,MOD(k+1+6,6),i,j),S2) 406 CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(k+1),:),S1) 407 CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(MOD(k+1+6,6)+1),:),S2) 143 408 Riv(n,k+1)=0.5*(S1+S2)/Ai(n) 144 409 Riv2(n,k+1)=0.5*(S1+S2)/surf_v(k+1) … … 236 501 237 502 ENDDO 503 CALL transfert_request(geom%Ai,req_i1) 504 CALL transfert_request(geom%centroid,req_i1) 238 505 CALL surf_triangle(d%xyz(:,ii_begin,jj_begin),d%xyz(:,ii_begin,jj_end),d%xyz(:,ii_end,jj_begin),S) 239 506 ! PRINT *,"Surf triangle : ",S*20/(4*Pi) -
codes/icosagcm/trunk/src/grid_param.f90
r12 r15 1 1 MODULE grid_param 2 INTEGER ,PARAMETER:: iim_glo=403 INTEGER ,PARAMETER :: jjm_glo=iim_glo2 INTEGER :: iim_glo=40 3 INTEGER :: jjm_glo 4 4 INTEGER,PARAMETER :: nb_face=10 5 INTEGER,PARAMETER :: llm=19 5 INTEGER :: llm=19 6 7 CONTAINS 8 9 SUBROUTINE init_grid_param 10 USE IOIPSL 11 IMPLICIT NONE 12 CALL getin('nbp',iim_glo) 13 jjm_glo=40 14 CALL getin('llm',llm) 15 16 END SUBROUTINE init_grid_param 17 6 18 END MODULE grid_param 7 19 -
codes/icosagcm/trunk/src/icosa_gcm.f90
r12 r15 9 9 USE timeloop_gcm_mod 10 10 USE transfert_mod 11 USE dissip_mod12 11 USE disvert_mod 13 12 USE etat0_mod 14 13 USE transfert_mod 14 USE vector 15 USE wind_mod 16 USE grid_param 15 17 IMPLICIT NONE 16 18 … … 20 22 INTEGER :: ind,i,j,k,n 21 23 REAL(rstd) :: tot_sum=0 24 REAL(rstd) :: vect(3,6) 25 REAL(rstd) :: centr(3),dist 22 26 27 CALL init_grid_param 23 28 CALL compute_metric 24 29 CALL compute_domain 30 CALL init_transfert 25 31 CALL compute_geometry 26 CALL init_transfert27 32 CALL init_disvert 28 33 … … 54 59 55 60 56 !CALL WriteField("Ai",geom%Ai)61 CALL WriteField("Ai",geom%Ai) 57 62 ! CALL WriteField("sum_ne",sum_ne) 58 59 63 60 64 CALL timeloop 61 65 -
codes/icosagcm/trunk/src/icosa_sw.f90
r13 r15 9 9 USE timeloop_sw_mod 10 10 USE transfert_mod 11 USE dissip_ mod11 USE dissip_sw_mod 12 12 USE disvert_mod 13 13 USE transfert_mod -
codes/icosagcm/trunk/src/metric.f90
r13 r15 128 128 ! PRINT *,"dist",vertex_glo(2,1,1)%xyz 129 129 ! PRINT *,"dist",p1/sqrt(sum(p1**2)) 130 CALL circ onscrit(vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,2,1)%xyz,p1)130 CALL circumcenter(vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,2,1)%xyz,p1) 131 131 ! CALL Centroide(vertex_glo(1,2,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,1,1)%xyz,p1) 132 132 ! CALL Centroide(vertex_glo(1,2,1)%xyz,vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,p1) -
codes/icosagcm/trunk/src/pression.f90
r12 r15 39 39 40 40 DO l = 1, llm+1 41 !$OMP DO 41 42 DO j=jj_begin-offset,jj_end+offset 42 43 DO i=ii_begin-offset,ii_end+offset -
codes/icosagcm/trunk/src/spherical_geom.f90
r13 r15 78 78 79 79 REAL(rstd) :: AB,AC,BC 80 REAL(rstd) :: s 80 REAL(rstd) :: s,x 81 81 82 82 CALL dist_cart(A,B,AB) … … 85 85 86 86 s=(AB+AC+BC)/2 87 surf=4*atan(sqrt( tan(s/2) * tan((s-AB)/2) * tan((s-AC)/2) * tan((s-BC)/2))) 87 x=tan(s/2) * tan((s-AB)/2) * tan((s-AC)/2) * tan((s-BC)/2) 88 IF (x<0) x=0. 89 surf=4*atan(sqrt( x)) 88 90 89 91 END SUBROUTINE surf_triangle … … 142 144 143 145 144 SUBROUTINE circ onscrit(A0,B0,C0,Center)146 SUBROUTINE circumcenter(A0,B0,C0,Center) 145 147 USE vector 146 148 IMPLICIT NONE … … 157 159 center=center/sqrt(sum(center**2)) 158 160 159 END SUBROUTINE circ onscrit161 END SUBROUTINE circumcenter 160 162 161 163 162 SUBROUTINE centroide(A,B,C,Center) 164 SUBROUTINE compute_centroid(points,n,centr) 165 USE vector 163 166 IMPLICIT NONE 164 REAL(rstd), INTENT(IN) :: A(3),B(3),C(3) 165 REAL(rstd), INTENT(OUT) :: Center(3) 167 INTEGER :: n 168 REAL(rstd), INTENT(IN) :: points(3,n) 169 REAL(rstd), INTENT(OUT) :: Centr(3) 170 171 REAL(rstd) :: p1(3),p2(3),cross(3) 172 REAL(rstd) :: norm_cross 173 INTEGER :: i,j 174 175 Centr(:)=0 176 DO i=1,n 177 j=MOD(i,n)+1 178 p1=points(:,i)/norm(points(:,i)) 179 p2=points(:,j)/norm(points(:,j)) 180 CALL cross_product2(p1,p2,cross) 181 norm_cross=norm(cross) 182 if (norm_cross<1e-10) CYCLE 183 184 Centr(:)=centr(:)+asin(norm_cross)*cross(:)/norm_cross 185 ENDDO 186 187 Centr(:)=centr(:)/norm(centr(:)) 166 188 167 REAL(rstd) :: AB05(3),AC05(3),BC05(3) 168 REAL(rstd) :: Vab(3),Vbc(3),Vca(3) 169 REAL(rstd) :: t 170 INTEGER :: n1,n2 171 172 AB05=(A+B)/2. 173 AC05=(A+C)/2. 174 BC05=(B+C)/2. 175 176 CALL director(A,B,C,Vab) 177 ! PRINT *,'director',sum((B-A)*Vab) 178 CALL director(B,C,A,Vbc) 179 ! PRINT *,'director',sum((C-B)*Vbc) 180 CALL director(C,A,B,Vca) 181 ! PRINT *,'director',sum((A-C)*Vca) 182 183 n1=1 ; n2=2 184 185 t=(Vbc(n1)*(BC05(n2)-AB05(n2))+(AB05(n1)-BC05(n1))*Vbc(n2))/(Vbc(n1)*Vab(n2)-Vab(n1)*Vbc(n2)) 186 Center=AB05+t*Vab 187 Center=Center/sqrt(sum(Center**2)) 188 ! PRINT*,"Center :",A 189 ! PRINT*,"Center :",B 190 ! PRINT*,"Center :",C 191 ! PRINT*,"Center :",Center 192 193 CONTAINS 194 195 SUBROUTINE director(A,B,C,V) 196 IMPLICIT NONE 197 REAL(rstd), INTENT(IN) :: A(3),B(3),C(3) 198 REAL(rstd), INTENT(OUT) :: V(3) 199 REAL(rstd) :: AB(3),AC(3) 200 REAL(rstd) :: t 201 202 AB=B-A 203 AC=C-A 204 205 t=sum(AB*AC)/sum(AB**2) 206 V=AB*t-AC 207 END SUBROUTINE director 208 209 END SUBROUTINE centroide 189 END SUBROUTINE compute_centroid 210 190 211 191 END MODULE spherical_geom_mod -
codes/icosagcm/trunk/src/theta_rhodz.f90
r12 r15 2 2 3 3 CONTAINS 4 5 SUBROUTINE theta_rhodz2theta(f_ps,f_theta_rhodz,f_theta) 6 USE transfert_mod 7 USE field_mod 8 USE dimensions 9 USE geometry 10 USE domain_mod 11 IMPLICIT NONE 12 TYPE(t_field), POINTER :: f_ps(:) 13 TYPE(t_field), POINTER :: f_theta_rhodz(:) 14 TYPE(t_field), POINTER :: f_theta(:) 15 16 REAL(rstd), POINTER :: ps(:) 17 REAL(rstd), POINTER :: theta_rhodz(:,:) 18 REAL(rstd), POINTER :: theta(:,:) 19 INTEGER :: ind 4 20 21 DO ind=1,ndomain 22 CALL swap_dimensions(ind) 23 CALL swap_geometry(ind) 24 ps=f_ps(ind) 25 theta_rhodz=f_theta_rhodz(ind) 26 theta=f_theta(ind) 27 CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0) 28 ENDDO 29 30 END SUBROUTINE theta_rhodz2theta 31 32 SUBROUTINE theta_rhodz2temperature(f_ps,f_theta_rhodz,f_temp) 33 USE transfert_mod 34 USE field_mod 35 USE dimensions 36 USE geometry 37 USE domain_mod 38 IMPLICIT NONE 39 TYPE(t_field), POINTER :: f_ps(:) 40 TYPE(t_field), POINTER :: f_theta_rhodz(:) 41 TYPE(t_field), POINTER :: f_temp(:) 42 43 REAL(rstd), POINTER :: ps(:) 44 REAL(rstd), POINTER :: theta_rhodz(:,:) 45 REAL(rstd), POINTER :: temp(:,:) 46 INTEGER :: ind 47 48 DO ind=1,ndomain 49 CALL swap_dimensions(ind) 50 CALL swap_geometry(ind) 51 ps=f_ps(ind) 52 theta_rhodz=f_theta_rhodz(ind) 53 temp=f_temp(ind) 54 CALL compute_theta_rhodz2temperature(ps, theta_rhodz,temp,0) 55 ENDDO 56 57 END SUBROUTINE theta_rhodz2temperature 58 5 59 SUBROUTINE theta2theta_rhodz(f_ps,f_theta,f_theta_rhodz) 6 60 USE transfert_mod … … 41 95 INTEGER,INTENT(IN) :: offset 42 96 INTEGER :: i,j,ij,l 43 REAL(rstd) :: p(iim*jjm,llm+1) 97 REAL(rstd),ALLOCATABLE,SAVE :: p(:,:) 98 99 !$OMP BARRIER 100 !$OMP MASTER 101 ALLOCATE( p(iim*jjm,llm+1)) 102 !$OMP END MASTER 103 !$OMP BARRIER 44 104 45 105 CALL compute_pression(ps,p,offset) 46 106 47 107 DO l = 1, llm 108 !$OMP DO 48 109 DO j=jj_begin-offset,jj_end+offset 49 110 DO i=ii_begin-offset,ii_end+offset … … 53 114 ENDDO 54 115 ENDDO 55 116 117 !$OMP BARRIER 118 !$OMP MASTER 119 DEALLOCATE( p) 120 !$OMP END MASTER 121 !$OMP BARRIER 122 56 123 END SUBROUTINE compute_theta2theta_rhodz 57 124 125 SUBROUTINE compute_theta_rhodz2theta(ps,theta_rhodz,theta,offset) 126 USE dimensions 127 USE geometry 128 USE metric 129 USE pression_mod 130 IMPLICIT NONE 131 REAL(rstd),INTENT(IN) :: ps(iim*jjm) 132 REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm) 133 REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm) 134 INTEGER,INTENT(IN) :: offset 135 INTEGER :: i,j,ij,l 136 REAL(rstd),SAVE,ALLOCATABLE :: p(:,:) 137 138 !$OMP BARRIER 139 !$OMP MASTER 140 ALLOCATE( p(iim*jjm,llm+1)) 141 !$OMP END MASTER 142 !$OMP BARRIER 143 CALL compute_pression(ps,p,offset) 144 145 DO l = 1, llm 146 DO j=jj_begin-offset,jj_end+offset 147 DO i=ii_begin-offset,ii_end+offset 148 ij=(j-1)*iim+i 149 theta(ij,l) = theta_rhodz(ij,l) / ((p(ij,l)-p(ij,l+1))/g) 150 ENDDO 151 ENDDO 152 ENDDO 153 154 !$OMP BARRIER 155 !$OMP MASTER 156 DEALLOCATE( p) 157 !$OMP END MASTER 158 !$OMP BARRIER 159 160 END SUBROUTINE compute_theta_rhodz2theta 161 162 SUBROUTINE compute_theta_rhodz2temperature(ps,theta_rhodz,temp,offset) 163 USE dimensions 164 USE geometry 165 USE metric 166 USE earth_const 167 USE pression_mod 168 USE exner_mod 169 IMPLICIT NONE 170 REAL(rstd),INTENT(IN) :: ps(iim*jjm) 171 REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm) 172 REAL(rstd),INTENT(OUT) :: temp(iim*jjm,llm) 173 INTEGER,INTENT(IN) :: offset 174 INTEGER :: i,j,ij,l 175 REAL(rstd) :: p(iim*jjm,llm+1) 176 REAL(rstd) :: pk(iim*jjm,llm) 177 REAL(rstd) :: pks(iim*jjm) 178 179 CALL compute_pression(ps,p,offset) 180 CALL compute_exner(ps,p,pks,pk,offset) 181 182 DO l = 1, llm 183 DO j=jj_begin-offset,jj_end+offset 184 DO i=ii_begin-offset,ii_end+offset 185 ij=(j-1)*iim+i 186 temp(ij,l) = ( theta_rhodz(ij,l) / ((p(ij,l)-p(ij,l+1))/g) ) * pk(ij,l) / cpp 187 ENDDO 188 ENDDO 189 ENDDO 190 191 END SUBROUTINE compute_theta_rhodz2temperature 192 58 193 END MODULE theta2theta_rhodz_mod -
codes/icosagcm/trunk/src/timeloop_gcm.f90
r12 r15 16 16 USE transfert_mod 17 17 USE metric 18 USE dissip_ mod18 USE dissip_gcm_mod 19 19 USE ioipsl 20 20 USE caldyn_gcm_mod … … 24 24 TYPE(t_field),POINTER :: f_phis(:) 25 25 TYPE(t_field),POINTER :: f_theta(:) 26 TYPE(t_field),POINTER :: f_dtheta(:) 26 27 TYPE(t_field),POINTER :: f_ps(:),f_psm1(:), f_psm2(:) 27 28 TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:) … … 58 59 CALL getin('matsuno_period',matsuno_period) 59 60 IF (scheme==leapfrog) matsuno_period=itaumax+1 60 61 62 61 63 62 CALL allocate_field(f_phis,field_t,type_real) … … 77 76 CALL allocate_field(f_dum2,field_u,type_real,llm) 78 77 78 CALL allocate_field(f_theta,field_t,type_real,llm) 79 CALL allocate_field(f_dtheta,field_t,type_real,llm) 80 79 81 CALL allocate_field(f_theta_rhodz,field_t,type_real,llm) 80 82 CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm) … … 84 86 CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm) 85 87 86 CALL allocate_caldyn 88 CALL init_dissip(dt) 89 CALL init_caldyn(dt) 87 90 88 91 ! CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u) … … 93 96 PRINT *,"It No :",It 94 97 95 CALL caldyn( f_phis,f_ps,f_theta_rhodz,f_u, f_dps, f_dtheta_rhodz, f_du)98 CALL caldyn(it,f_phis,f_ps,f_theta_rhodz,f_u, f_dps, f_dtheta_rhodz, f_du) 96 99 97 100 IF (scheme==Euler) THEN … … 102 105 CALL leapfrog_matsuno_scheme 103 106 ELSE IF (scheme==adam_bashforth) THEN 107 CALL dissip(f_u,f_du,f_ps,f_theta_rhodz,f_dtheta_rhodz) 104 108 CALL adam_bashforth_scheme 105 109 ENDIF -
codes/icosagcm/trunk/src/timeloop_sw.f90
r13 r15 18 18 USE transfert_mod 19 19 USE metric 20 USE dissip_ mod20 USE dissip_sw_mod 21 21 USE ioipsl 22 22 … … 87 87 CALL leapfrog_matsuno_scheme 88 88 ELSE IF (scheme==adam_bashforth) THEN 89 CALL dissip(f_u,f_du) 89 90 CALL adam_bashforth_scheme 90 91 ENDIF 91 92 CALL dissip(f_u)93 92 94 93 ENDDO -
codes/icosagcm/trunk/src/transfert.f90
r12 r15 31 31 32 32 DO ind=1,ndomain 33 CALL swap_dimensions(ind) 33 34 DO i=ii_begin,ii_end+1 34 35 CALL request_add_point(ind,i,jj_begin-1,req_i1) … … 67 68 CALL create_request(field_u,req_e1) 68 69 DO ind=1,ndomain 70 CALL swap_dimensions(ind) 69 71 DO i=ii_begin,ii_end 70 72 CALL request_add_point(ind,i,jj_begin-1,req_e1,rup) -
codes/icosagcm/trunk/src/write_field.f90
r12 r15 402 402 n=n+1 403 403 CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n)) 404 lon(n)=lon(n)*180/Pi +180404 lon(n)=lon(n)*180/Pi 405 405 lat(n)=lat(n)*180/Pi 406 406 DO k=0,5 407 407 CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n)) 408 408 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 409 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi +180409 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 410 410 ENDDO 411 411 ENDIF … … 505 505 n=n+1 506 506 CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n)) 507 lon(n)=lon(n)*180/Pi+180 507 lon(n)=lon(n)*180/Pi 508 ! IF (lon(n)<0) lon(n)=lon(n)+360 508 509 lat(n)=lat(n)*180/Pi 509 510 CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) … … 513 514 DO k=0,2 514 515 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 515 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi+180 516 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 517 ! IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360 516 518 ENDDO 517 519 ENDDO … … 523 525 n=n+1 524 526 CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n)) 525 lon(n)=lon(n)*180/Pi+180 527 lon(n)=lon(n)*180/Pi 528 ! IF (lon(n)<0) lon(n)=lon(n)+360 526 529 lat(n)=lat(n)*180/Pi 527 530 CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) … … 531 534 DO k=0,2 532 535 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 533 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi+180 536 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 537 ! IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360 534 538 ENDDO 535 539 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.