![(please configure the [header_logo] section in trac.ini)](/dynamico/chrome/site/your_project_logo.png){kind=logo}
Changeset 267 for codes/icosagcm/branches
- Timestamp:
- 08/05/14 16:00:09 (10 years ago)
- Location:
- codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src
- Files:
-
- 9 edited
-
advect.f90 (modified) (6 diffs)
-
advect_tracer.f90 (modified) (6 diffs)
-
caldyn.f90 (modified) (2 diffs)
-
dissip_gcm.f90 (modified) (1 diff)
-
getin.f90 (modified) (8 diffs)
-
mpipara.F90 (modified) (1 diff)
-
spherical_geom.f90 (modified) (4 diffs)
-
time.f90 (modified) (2 diffs)
-
transfert_omp.f90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/advect.f90
r221 r267 11 11 !========================================================================== 12 12 13 SUBROUTINE init_advect(normal,tangent, one_over_sqrt_leng)13 SUBROUTINE init_advect(normal,tangent,sqrt_leng) 14 14 IMPLICIT NONE 15 15 REAL(rstd),INTENT(OUT) :: normal(3*iim*jjm,3) 16 16 REAL(rstd),INTENT(OUT) :: tangent(3*iim*jjm,3) 17 REAL(rstd),INTENT(OUT) :: one_over_sqrt_leng(iim*jjm)17 REAL(rstd),INTENT(OUT) :: sqrt_leng(iim*jjm) 18 18 19 19 INTEGER :: ij … … 40 40 tangent(ij+u_ldown,:)=tangent(ij+u_ldown,:)/sqrt(sum(tangent(ij+u_ldown,:)**2)+1e-50) 41 41 42 one_over_sqrt_leng(ij) = 1./sqrt(max(sum((xyz_v(ij+z_up,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_down,:) - xyz_i(ij,:))**2), &43 44 42 sqrt_leng(ij) = sqrt(max(sum((xyz_v(ij+z_up,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_down,:) - xyz_i(ij,:))**2), & 43 sum((xyz_v(ij+z_rup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_rdown,:) - xyz_i(ij,:))**2), & 44 sum((xyz_v(ij+z_lup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_ldown,:) - xyz_i(ij,:))**2)) ) 45 45 ENDDO 46 46 … … 49 49 !======================================================================================= 50 50 51 SUBROUTINE compute_gradq3d(qi_in, one_over_sqrt_leng_in,gradq3d_out,xyz_i,xyz_v)51 SUBROUTINE compute_gradq3d(qi_in,sqrt_leng_in,gradq3d_out,xyz_i,xyz_v) 52 52 USE trace 53 53 USE omp_para 54 54 IMPLICIT NONE 55 55 REAL(rstd),INTENT(IN) :: qi_in(iim*jjm,llm) 56 REAL(rstd),INTENT(IN) :: one_over_sqrt_leng_in(iim*jjm)56 REAL(rstd),INTENT(IN) :: sqrt_leng_in(iim*jjm) 57 57 REAL(rstd),INTENT(IN) :: xyz_i(iim*jjm,3) 58 58 REAL(rstd),INTENT(IN) :: xyz_v(2*iim*jjm,3) … … 66 66 INTEGER :: ij,k,ind,l 67 67 REAL(rstd) :: qi(iim*jjm,llm) 68 REAL(rstd) :: one_over_sqrt_leng(iim*jjm)68 REAL(rstd) :: sqrt_leng(iim*jjm) 69 69 REAL(rstd) :: gradq3d(iim*jjm,llm,3) 70 70 REAL(rstd) :: detx,dety,detz,det … … 74 74 75 75 qi=qi_in 76 one_over_sqrt_leng=one_over_sqrt_leng_in76 sqrt_leng=sqrt_leng_in 77 77 78 78 CALL trace_start("compute_gradq3d1") … … 252 252 maggrd = gradq3d(ij,l,1)*gradq3d(ij,l,1) + gradq3d(ij,l,2)*gradq3d(ij,l,2) + gradq3d(ij,l,3)*gradq3d(ij,l,3) 253 253 maggrd = sqrt(maggrd) 254 maxq_c = qi(ij,l) + maggrd* one_over_sqrt_leng(ij)255 minq_c = qi(ij,l) - maggrd* one_over_sqrt_leng(ij)254 maxq_c = qi(ij,l) + maggrd*sqrt_leng(ij) 255 minq_c = qi(ij,l) - maggrd*sqrt_leng(ij) 256 256 maxq = max(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), & 257 257 qi(ij+t_rdown,l),qi(ij+t_ldown,l)) -
codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/advect_tracer.f90
r221 r267 8 8 TYPE(t_field),SAVE,POINTER :: f_gradq3d(:) 9 9 TYPE(t_field),SAVE,POINTER :: f_cc(:) ! starting point of backward-trajectory (Miura approach) 10 TYPE(t_field),SAVE,POINTER :: f_ one_over_sqrt_leng(:)10 TYPE(t_field),SAVE,POINTER :: f_sqrt_leng(:) 11 11 12 12 TYPE(t_message),SAVE :: req_u, req_cc, req_wfluxt, req_q, req_rhodz, req_gradq3d … … 28 28 REAL(rstd),POINTER :: tangent(:,:) 29 29 REAL(rstd),POINTER :: normal(:,:) 30 REAL(rstd),POINTER :: one_over_sqrt_leng(:)30 REAL(rstd),POINTER :: sqrt_leng(:) 31 31 INTEGER :: ind 32 32 … … 35 35 CALL allocate_field(f_gradq3d,field_t,type_real,llm,3, name='gradq3d') 36 36 CALL allocate_field(f_cc,field_u,type_real,llm,3, name='cc') 37 CALL allocate_field(f_ one_over_sqrt_leng,field_t,type_real, name='one_over_sqrt_leng')37 CALL allocate_field(f_sqrt_leng,field_t,type_real, name='sqrt_leng') 38 38 CALL allocate_field(f_dzqw, field_t, type_real, llm, name='dzqw') 39 39 CALL allocate_field(f_adzqw, field_t, type_real, llm, name='adzqw') … … 47 47 normal=f_normal(ind) 48 48 tangent=f_tangent(ind) 49 one_over_sqrt_leng=f_one_over_sqrt_leng(ind)50 CALL init_advect(normal,tangent, one_over_sqrt_leng)49 sqrt_leng=f_sqrt_leng(ind) 50 CALL init_advect(normal,tangent,sqrt_leng) 51 51 END DO 52 52 … … 66 66 TYPE(t_field),POINTER :: f_rhodz(:) ! mass field at beginning of macro time step 67 67 68 REAL(rstd),POINTER :: q(:,:,:), normal(:,:), tangent(:,:), one_over_sqrt_leng(:), gradq3d(:,:,:), cc(:,:,:)68 REAL(rstd),POINTER :: q(:,:,:), normal(:,:), tangent(:,:), sqrt_leng(:), gradq3d(:,:,:), cc(:,:,:) 69 69 REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:) 70 70 REAL(rstd),POINTER :: rhodz(:,:), u(:,:) … … 145 145 q = f_q(ind) 146 146 gradq3d = f_gradq3d(ind) 147 one_over_sqrt_leng=f_one_over_sqrt_leng(ind)148 CALL compute_gradq3d(q(:,:,k), one_over_sqrt_leng,gradq3d,xyz_i,xyz_v)147 sqrt_leng=f_sqrt_leng(ind) 148 CALL compute_gradq3d(q(:,:,k),sqrt_leng,gradq3d,xyz_i,xyz_v) 149 149 END DO 150 150 -
codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/caldyn.f90
r221 r267 1 1 MODULE caldyn_mod 2 2 USE icosa 3 USE caldyn_gcm_mod, ONLY : caldyn_BC4 3 PRIVATE 5 4 CHARACTER(LEN=255),SAVE :: caldyn_type … … 31 30 32 31 END SUBROUTINE init_caldyn 32 33 SUBROUTINE caldyn_BC(f_phis, f_wflux) 34 USE caldyn_gcm_mod, ONLY : caldyn_gcm_BC=>caldyn_BC 35 IMPLICIT NONE 36 TYPE(t_field), POINTER :: f_phis(:), f_wflux(:) 37 SELECT CASE (TRIM(caldyn_type)) 38 CASE('gcm') 39 CALL caldyn_gcm_BC(f_phis, f_wflux) 40 END SELECT 41 END SUBROUTINE caldyn_BC 33 42 34 43 SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, & -
codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/dissip_gcm.f90
r221 r267 427 427 ENDDO 428 428 429 itau_dissip=INT(mintau/dt) 429 IF(mintau>0) THEN 430 itau_dissip=INT(mintau/dt) 431 dtdissip=itau_dissip*dt 432 ELSE 433 IF (is_mpi_root) PRINT *,"No dissipation time set, setting itau_dissip to 1000000000" 434 itau_dissip=100000000 435 END IF 430 436 itau_dissip=MAX(1,itau_dissip) 431 dtdissip=itau_dissip*dt432 437 IF (is_mpi_root) PRINT *,"mintau ",mintau,"itau_dissip",itau_dissip," dtdissip ",dtdissip 433 438 -
codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/getin.f90
r221 r267 15 15 USE transfert_omp_mod 16 16 USE omp_para 17 USE mpipara 17 18 IMPLICIT NONE 18 19 CHARACTER(LEN=*) :: name … … 21 22 !$OMP MASTER 22 23 CALL getin_(name,value) 24 IF(is_mpi_root) PRINT *,'GETIN ',TRIM(name),' = ', TRIM(value) 23 25 !$OMP END MASTER 24 26 IF (omp_in_parallel()) CALL bcast_omp(value) … … 29 31 USE transfert_omp_mod 30 32 USE omp_para 33 USE mpipara 31 34 IMPLICIT NONE 32 35 CHARACTER(LEN=*) :: name … … 35 38 !$OMP MASTER 36 39 CALL getin_(name,value) 40 IF(is_mpi_root) PRINT *,'GETIN ',TRIM(name),' = ', value 37 41 !$OMP END MASTER 38 42 IF (omp_in_parallel()) CALL bcast_omp(value) … … 43 47 USE transfert_omp_mod 44 48 USE omp_para 49 USE mpipara 45 50 IMPLICIT NONE 46 51 CHARACTER(LEN=*) :: name … … 49 54 !$OMP MASTER 50 55 CALL getin_(name,value) 56 IF(is_mpi_root) PRINT *,'GETIN ',TRIM(name), ' = ', value 51 57 !$OMP END MASTER 52 58 IF (omp_in_parallel()) CALL bcast_omp(value) … … 57 63 USE ioipsl, ONLY : getin_=>getin 58 64 USE omp_para 65 USE mpipara 59 66 USE transfert_omp_mod 60 67 IMPLICIT NONE … … 64 71 !$OMP MASTER 65 72 CALL getin_(name,value) 73 IF(is_mpi_root) PRINT *,'GETIN ',TRIM(name), ' = ', value 66 74 !$OMP END MASTER 67 75 IF (omp_in_parallel()) CALL bcast_omp(value) -
codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/mpipara.F90
r261 r267 21 21 END INTERFACE free_mpi_buffer 22 22 23 PRIVATE :: getin 24 23 25 CONTAINS 24 26 27 SUBROUTINE getin(name,value) ! Copied from getin.f90 to avoid circular dependency 28 USE ioipsl, ONLY : getin_=>getin 29 USE transfert_omp_mod 30 USE omp_para 31 IMPLICIT NONE 32 CHARACTER(LEN=*) :: name 33 CHARACTER(LEN=*) :: value 34 35 !$OMP MASTER 36 CALL getin_(name,value) 37 IF(is_mpi_root) PRINT *,'GETIN ',TRIM(name),' = ', TRIM(value) 38 !$OMP END MASTER 39 IF (omp_in_parallel()) CALL bcast_omp(value) 40 END SUBROUTINE getin 41 25 42 SUBROUTINE init_mpipara 26 43 USE mpi_mod 27 USE getin_mod28 44 #ifdef CPP_USING_XIOS 29 45 USE xios -
codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/spherical_geom.f90
r221 r267 179 179 180 180 181 SUBROUTINE circumcenter(A0,B0,C0, Center)181 SUBROUTINE circumcenter(A0,B0,C0,center) 182 182 USE vector 183 183 IMPLICIT NONE … … 185 185 REAL(rstd), INTENT(OUT) :: Center(3) 186 186 187 REAL(rstd) :: a(3),b(3),c(3) 187 REAL(rstd) :: a(3),b(3),c(3), ac(3), ab(3), p1(3), q(3), p2(3) 188 188 189 189 a=A0/sqrt(sum(A0**2)) … … 191 191 c=C0/sqrt(sum(C0**2)) 192 192 193 CALL Cross_product2(b-a,c-b,center) 194 center=center/sqrt(sum(center**2)) 195 193 ab=b-a 194 ac=c-a 195 CALL Cross_product2(ab,ac,p1) 196 IF(.FALSE.) THEN ! Direct solution, round-off error 197 center=p1/norm(p1) 198 ELSE ! Two-step solution, stable 199 q = SUM(ac**2)*ab-SUM(ab**2)*ac 200 CALL Cross_product2(p1,q,p2) 201 p2 = a + p2/(2.*SUM(p1**2)) 202 center = p2/norm(p2) 203 END IF 196 204 END SUBROUTINE circumcenter 197 205 … … 204 212 REAL(rstd), INTENT(OUT) :: Centr(3) 205 213 206 REAL(rstd) :: p1(3),p2(3),cross(3) 207 REAL(rstd) :: norm_cross 214 REAL(rstd) :: p1(3),p2(3),cross(3), cc(3) 215 REAL(rstd) :: norm_cross, area 208 216 INTEGER :: i,j 209 210 Centr(:)=0 211 DO i=1,n 212 j=MOD(i,n)+1 213 p1=points(:,i)/norm(points(:,i)) 214 p2=points(:,j)/norm(points(:,j)) 215 CALL cross_product2(p1,p2,cross) 216 norm_cross=norm(cross) 217 if (norm_cross<1e-10) CYCLE 218 219 Centr(:)=centr(:)+asin(norm_cross)*cross(:)/norm_cross 220 ENDDO 221 222 Centr(:)=centr(:)/norm(centr(:)) 223 217 218 centr(:)=0 219 IF(.FALSE.) THEN 220 ! Gauss formula (subject to round-off error) 221 DO i=1,n 222 j=MOD(i,n)+1 223 p1=points(:,i)/norm(points(:,i)) 224 p2=points(:,j)/norm(points(:,j)) 225 CALL cross_product2(p1,p2,cross) 226 norm_cross=norm(cross) 227 if (norm_cross<1e-10) CYCLE 228 centr(:)=centr(:)+asin(norm_cross)*cross(:)/norm_cross 229 ENDDO 230 ELSE 231 ! Simple area-weighted average (second-order accurate, stable) 232 cc=SUM(points,2) ! arithmetic average used as center 233 cc=cc/norm(cc) 234 DO i=1,n 235 j=MOD(i,n)+1 236 p1=points(:,i)/norm(points(:,i)) 237 p2=points(:,j)/norm(points(:,j)) 238 CALL surf_triangle(cc,p1,p2,area) 239 centr(:)=centr(:)+area*(p1+p2+cc) 240 ENDDO 241 END IF 242 243 centr(:)=centr(:)/norm(centr(:)) 244 224 245 END SUBROUTINE compute_centroid 225 246 -
codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/time.f90
r262 r267 68 68 CALL getin('run_length',run_length) 69 69 itaumax=run_length/dt 70 IF (is_mpi_root) THEN71 PRINT *, 'itaumax=',itaumax72 PRINT *, 'itau_adv=',itau_adv, 'itau_dissip=',itau_dissip, 'itau_physics=',itau_physics73 END IF74 70 dt=dt/scale_factor 75 71 … … 90 86 CALL getin('itau_physics',itau_physics) 91 87 88 IF (is_mpi_root) THEN 89 PRINT *, 'itaumax=',itaumax 90 PRINT *, 'itau_adv=',itau_adv, 'itau_dissip=',itau_dissip, 'itau_physics=',itau_physics 91 END IF 92 92 93 93 CALL create_time_counter_header -
codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/transfert_omp.f90
r221 r267 95 95 IMPLICIT NONE 96 96 CHARACTER(LEN=*),INTENT(INOUT) :: Var 97 98 CALL bcast_omp_cgen(Var,len(Var),buffer_c) 99 97 INTEGER :: lenvar 98 lenvar=len(Var) 99 CALL bcast_omp_i(lenvar) 100 CALL bcast_omp_cgen(Var,lenvar,buffer_c) 100 101 END SUBROUTINE bcast_omp_c 101 102
Note: See TracChangeset
for help on using the changeset viewer.
