Changeset 186 for codes/icosagcm/trunk
- Timestamp:
- 01/09/14 09:56:11 (10 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 2 added
- 58 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/advect.f90
r174 r186 49 49 !======================================================================================= 50 50 51 SUBROUTINE compute_gradq3d(qi ,one_over_sqrt_leng,gradq3d)51 SUBROUTINE compute_gradq3d(qi_in,one_over_sqrt_leng_in,gradq3d_out,xyz_i,xyz_v) 52 52 USE trace 53 53 USE omp_para 54 54 IMPLICIT NONE 55 REAL(rstd),INTENT(IN) :: qi(iim*jjm,llm) 56 REAL(rstd),INTENT(IN) :: one_over_sqrt_leng(iim*jjm) 57 REAL(rstd),INTENT(OUT) :: gradq3d(iim*jjm,llm,3) 55 REAL(rstd),INTENT(IN) :: qi_in(iim*jjm,llm) 56 REAL(rstd),INTENT(IN) :: one_over_sqrt_leng_in(iim*jjm) 57 REAL(rstd),INTENT(IN) :: xyz_i(iim*jjm,3) 58 REAL(rstd),INTENT(IN) :: xyz_v(2*iim*jjm,3) 59 REAL(rstd),INTENT(OUT) :: gradq3d_out(iim*jjm,llm,3) 58 60 REAL(rstd) :: maxq,minq,minq_c,maxq_c 59 61 REAL(rstd) :: alphamx,alphami,alpha ,maggrd … … 63 65 REAL(rstd) :: gradtri(2*iim*jjm,llm,3) 64 66 INTEGER :: ij,k,ind,l 65 66 CALL trace_start("compute_gradq3d") 67 REAL(rstd) :: qi(iim*jjm,llm) 68 REAL(rstd) :: one_over_sqrt_leng(iim*jjm) 69 REAL(rstd) :: gradq3d(iim*jjm,llm,3) 70 REAL(rstd) :: detx,dety,detz,det 71 REAL(rstd) :: A(3,3), a11,a12,a13,a21,a22,a23,a31,a32,a33 72 REAL(rstd) :: x1,x2,x3 73 REAL(rstd) :: dq(3) 74 75 qi=qi_in 76 one_over_sqrt_leng=one_over_sqrt_leng_in 77 78 CALL trace_start("compute_gradq3d1") 67 79 68 80 ! TODO : precompute ar, drop arr as output argument of gradq ? … … 71 83 ! Compute gradient at triangles solving a linear system 72 84 ! arr = area of triangle joining centroids of hexagons 85 ! DO l = ll_begin,ll_end 86 !!$SIMD 87 ! DO ij=ij_begin_ext,ij_end_ext 88 !! CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,:),arr(ij+z_up)) 89 !! CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,:),arr(ij+z_down)) 90 ! CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,1),gradtri(ij+z_up,l,2),gradtri(ij+z_up,l,3),arr(ij+z_up)) 91 ! CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,1),gradtri(ij+z_down,l,2),gradtri(ij+z_down,l,3),arr(ij+z_down)) 92 ! END DO 93 ! END DO 94 73 95 DO l = ll_begin,ll_end 74 96 !$SIMD 75 97 DO ij=ij_begin_ext,ij_end_ext 76 CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,:),arr(ij+z_up)) 77 CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,:),arr(ij+z_down)) 78 END DO 98 ! CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,1),gradtri(ij+z_up,l,2),gradtri(ij+z_up,l,3),arr(ij+z_up)) 99 100 101 A(1,1)=xyz_i(ij+t_rup,1)-xyz_i(ij,1); A(1,2)=xyz_i(ij+t_rup,2)-xyz_i(ij,2); A(1,3)=xyz_i(ij+t_rup,3)-xyz_i(ij,3) 102 A(2,1)=xyz_i(ij+t_lup,1)-xyz_i(ij,1); A(2,2)=xyz_i(ij+t_lup,2)-xyz_i(ij,2); A(2,3)=xyz_i(ij+t_lup,3)-xyz_i(ij,3) 103 A(3,1)=xyz_v(ij+z_up,1); A(3,2)= xyz_v(ij+z_up,2); A(3,3)=xyz_v(ij+z_up,3) 104 105 dq(1) = qi(ij+t_rup,l)-qi(ij,l) 106 dq(2) = qi(ij+t_lup,l)-qi(ij,l) 107 dq(3) = 0.0 108 109 110 ! CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det) 111 112 a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1) 113 a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2) 114 a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3) 115 116 x1 = a11 * (a22 * a33 - a23 * a32) 117 x2 = a12 * (a21 * a33 - a23 * a31) 118 x3 = a13 * (a21 * a32 - a22 * a31) 119 det = x1 - x2 + x3 120 121 ! CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),detx) 122 123 a11=dq(1) ; a12=dq(2) ; a13=dq(3) 124 a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2) 125 a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3) 126 127 x1 = a11 * (a22 * a33 - a23 * a32) 128 x2 = a12 * (a21 * a33 - a23 * a31) 129 x3 = a13 * (a21 * a32 - a22 * a31) 130 detx = x1 - x2 + x3 131 132 ! CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dety) 133 134 a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1) 135 a21=dq(1) ; a22=dq(2) ; a23=dq(3) 136 a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3) 137 138 x1 = a11 * (a22 * a33 - a23 * a32) 139 x2 = a12 * (a21 * a33 - a23 * a31) 140 x3 = a13 * (a21 * a32 - a22 * a31) 141 dety = x1 - x2 + x3 142 143 ! CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),detz) 144 145 a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1) 146 a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2) 147 a31=dq(1) ; a32=dq(2) ; a33=dq(3) 148 149 x1 = a11 * (a22 * a33 - a23 * a32) 150 x2 = a12 * (a21 * a33 - a23 * a31) 151 x3 = a13 * (a21 * a32 - a22 * a31) 152 detz = x1 - x2 + x3 153 154 gradtri(ij+z_up,l,1) = detx 155 gradtri(ij+z_up,l,2) = dety 156 gradtri(ij+z_up,l,3) = detz 157 arr(ij+z_up) = det 158 159 ENDDO 160 161 DO ij=ij_begin_ext,ij_end_ext 162 163 164 ! CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,1),gradtri(ij+z_down,l,2),gradtri(ij+z_down,l,3),arr(ij+z_down)) 165 166 A(1,1)=xyz_i(ij+t_ldown,1)-xyz_i(ij,1); A(1,2)=xyz_i(ij+t_ldown,2)-xyz_i(ij,2); A(1,3)=xyz_i(ij+t_ldown,3)-xyz_i(ij,3) 167 A(2,1)=xyz_i(ij+t_rdown,1)-xyz_i(ij,1); A(2,2)=xyz_i(ij+t_rdown,2)-xyz_i(ij,2); A(2,3)=xyz_i(ij+t_rdown,3)-xyz_i(ij,3) 168 A(3,1)=xyz_v(ij+z_down,1); A(3,2)= xyz_v(ij+z_down,2); A(3,3)=xyz_v(ij+z_down,3) 169 170 dq(1) = qi(ij+t_ldown,l)-qi(ij,l) 171 dq(2) = qi(ij+t_rdown,l)-qi(ij,l) 172 dq(3) = 0.0 173 174 175 ! CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det) 176 177 a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1) 178 a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2) 179 a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3) 180 181 x1 = a11 * (a22 * a33 - a23 * a32) 182 x2 = a12 * (a21 * a33 - a23 * a31) 183 x3 = a13 * (a21 * a32 - a22 * a31) 184 det = x1 - x2 + x3 185 186 ! CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),detx) 187 188 a11=dq(1) ; a12=dq(2) ; a13=dq(3) 189 a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2) 190 a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3) 191 192 x1 = a11 * (a22 * a33 - a23 * a32) 193 x2 = a12 * (a21 * a33 - a23 * a31) 194 x3 = a13 * (a21 * a32 - a22 * a31) 195 detx = x1 - x2 + x3 196 197 ! CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dety) 198 199 a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1) 200 a21=dq(1) ; a22=dq(2) ; a23=dq(3) 201 a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3) 202 203 x1 = a11 * (a22 * a33 - a23 * a32) 204 x2 = a12 * (a21 * a33 - a23 * a31) 205 x3 = a13 * (a21 * a32 - a22 * a31) 206 dety = x1 - x2 + x3 207 208 ! CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),detz) 209 210 a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1) 211 a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2) 212 a31=dq(1) ; a32=dq(2) ; a33=dq(3) 213 214 x1 = a11 * (a22 * a33 - a23 * a32) 215 x2 = a12 * (a21 * a33 - a23 * a31) 216 x3 = a13 * (a21 * a32 - a22 * a31) 217 detz = x1 - x2 + x3 218 219 gradtri(ij+z_down,l,1) = detx 220 gradtri(ij+z_down,l,2) = dety 221 gradtri(ij+z_down,l,3) = detz 222 arr(ij+z_down) = det 223 224 END DO 79 225 END DO 80 226 81 227 !$SIMD 82 83 228 DO ij=ij_begin,ij_end 229 ar(ij) = arr(ij+z_up)+arr(ij+z_lup)+arr(ij+z_ldown)+arr(ij+z_down)+arr(ij+z_rdown)+arr(ij+z_rup)+1.e-50 84 230 ENDDO 231 CALL trace_end("compute_gradq3d1") 232 CALL trace_start2("compute_gradq3d2") 85 233 86 234 DO k=1,3 … … 94 242 END DO 95 243 ENDDO 244 CALL trace_end2("compute_gradq3d2") 245 CALL trace_start("compute_gradq3d3") 96 246 97 247 !============================================================================================= LIMITING … … 99 249 !$SIMD 100 250 DO ij=ij_begin,ij_end 101 maggrd = dot_product(gradq3d(ij,l,:),gradq3d(ij,l,:)) 251 ! maggrd = dot_product(gradq3d(ij,l,:),gradq3d(ij,l,:)) 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) 102 253 maggrd = sqrt(maggrd) 103 254 maxq_c = qi(ij,l) + maggrd*one_over_sqrt_leng(ij) … … 112 263 alphami = max(alphami,0.0) 113 264 alpha = min(alphamx,alphami,1.0) 114 gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:) 265 ! gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:) 266 gradq3d(ij,l,1) = alpha*gradq3d(ij,l,1) 267 gradq3d(ij,l,2) = alpha*gradq3d(ij,l,2) 268 gradq3d(ij,l,3) = alpha*gradq3d(ij,l,3) 115 269 END DO 116 270 END DO 117 271 118 CALL trace_end("compute_gradq3d") 272 CALL trace_end("compute_gradq3d3") 273 274 gradq3d_out=gradq3d 275 276 CONTAINS 277 278 SUBROUTINE gradq(n0,l,n1,n2,n3,q,dq1,dq2,dq3,det) 279 IMPLICIT NONE 280 INTEGER, INTENT(IN) :: n0,l,n1,n2,n3 281 REAL(rstd), INTENT(IN) :: q(iim*jjm,llm) 282 ! REAL(rstd), INTENT(OUT) :: dq(3), det 283 REAL(rstd), INTENT(OUT) :: dq1,dq2,dq3,det 284 REAL(rstd) :: dq(3) 285 286 REAL(rstd) ::detx,dety,detz 287 INTEGER :: info 288 INTEGER :: IPIV(3) 289 290 REAL(rstd) :: A(3,3) 291 REAL(rstd) :: B(3) 292 293 ! TODO : replace A by A1,A2,A3 294 A(1,1)=xyz_i(n1,1)-xyz_i(n0,1); A(1,2)=xyz_i(n1,2)-xyz_i(n0,2); A(1,3)=xyz_i(n1,3)-xyz_i(n0,3) 295 A(2,1)=xyz_i(n2,1)-xyz_i(n0,1); A(2,2)=xyz_i(n2,2)-xyz_i(n0,2); A(2,3)=xyz_i(n2,3)-xyz_i(n0,3) 296 A(3,1)=xyz_v(n3,1); A(3,2)= xyz_v(n3,2); A(3,3)=xyz_v(n3,3) 297 298 dq(1) = q(n1,l)-q(n0,l) 299 dq(2) = q(n2,l)-q(n0,l) 300 dq(3) = 0.0 301 302 ! CALL DGESV(3,1,A,3,IPIV,dq(:),3,info) 303 ! CALL determinant(A(:,1),A(:,2),A(:,3),det) 304 ! CALL determinant(dq,A(:,2),A(:,3),detx) 305 ! CALL determinant(A(:,1),dq,A(:,3),dety) 306 ! CALL determinant(A(:,1),A(:,2),dq,detz) 307 ! dq(1) = detx 308 ! dq(2) = dety 309 ! dq(3) = detz 310 311 CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det) 312 CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),dq1) 313 CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dq2) 314 CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),dq3) 315 316 END SUBROUTINE gradq 317 318 !========================================================================== 319 ! PURE SUBROUTINE determinant(a1,a2,a3,det) 320 ! IMPLICIT NONE 321 ! REAL(rstd), DIMENSION(3), INTENT(IN) :: a1,a2,a3 322 ! REAL(rstd), INTENT(OUT) :: det 323 ! REAL(rstd) :: x1,x2,x3 324 ! x1 = a1(1) * (a2(2) * a3(3) - a2(3) * a3(2)) 325 ! x2 = a1(2) * (a2(1) * a3(3) - a2(3) * a3(1)) 326 ! x3 = a1(3) * (a2(1) * a3(2) - a2(2) * a3(1)) 327 ! det = x1 - x2 + x3 328 ! END SUBROUTINE determinant 329 330 SUBROUTINE determinant(a11,a12,a13,a21,a22,a23,a31,a32,a33,det) 331 IMPLICIT NONE 332 REAL(rstd), INTENT(IN) :: a11,a12,a13,a21,a22,a23,a31,a32,a33 333 REAL(rstd), INTENT(OUT) :: det 334 REAL(rstd) :: x1,x2,x3 335 x1 = a11 * (a22 * a33 - a23 * a32) 336 x2 = a12 * (a21 * a33 - a23 * a31) 337 x3 = a13 * (a21 * a32 - a22 * a31) 338 det = x1 - x2 + x3 339 END SUBROUTINE determinant 340 341 119 342 END SUBROUTINE compute_gradq3d 120 343 … … 222 445 IF (ne(ij,right)*hfluxt(ij+u_right,l)>0) THEN 223 446 ed = cc(ij+u_right,l,:) - xyz_i(ij,:) 224 qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 447 ! qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 448 qe = qi(ij,l)+gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 225 449 ELSE 226 450 ed = cc(ij+u_right,l,:) - xyz_i(ij+t_right,:) 227 qe = qi(ij+t_right,l)+sum2(gradq3d(ij+t_right,l,:),ed) 451 ! qe = qi(ij+t_right,l)+sum2(gradq3d(ij+t_right,l,:),ed) 452 qe = qi(ij+t_right,l) + gradq3d(ij+t_right,l,1)*ed(1)+gradq3d(ij+t_right,l,2)*ed(2)+gradq3d(ij+t_right,l,3)*ed(3) 228 453 ENDIF 229 454 qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe … … 231 456 IF (ne(ij,lup)*hfluxt(ij+u_lup,l)>0) THEN 232 457 ed = cc(ij+u_lup,l,:) - xyz_i(ij,:) 233 qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 458 ! qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 459 qe = qi(ij,l) + gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 234 460 ELSE 235 461 ed = cc(ij+u_lup,l,:) - xyz_i(ij+t_lup,:) 236 qe = qi(ij+t_lup,l)+sum2(gradq3d(ij+t_lup,l,:),ed) 462 ! qe = qi(ij+t_lup,l)+sum2(gradq3d(ij+t_lup,l,:),ed) 463 qe = qi(ij+t_lup,l) + gradq3d(ij+t_lup,l,1)*ed(1)+gradq3d(ij+t_lup,l,2)*ed(2)+gradq3d(ij+t_lup,l,3)*ed(3) 237 464 ENDIF 238 465 qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe … … 240 467 IF (ne(ij,ldown)*hfluxt(ij+u_ldown,l)>0) THEN 241 468 ed = cc(ij+u_ldown,l,:) - xyz_i(ij,:) 242 qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 469 ! qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 470 qe = qi(ij,l) + gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 243 471 ELSE 244 472 ed = cc(ij+u_ldown,l,:) - xyz_i(ij+t_ldown,:) 245 qe = qi(ij+t_ldown,l)+sum2(gradq3d(ij+t_ldown,l,:),ed) 473 ! qe = qi(ij+t_ldown,l)+sum2(gradq3d(ij+t_ldown,l,:),ed) 474 qe = qi(ij+t_ldown,l)+gradq3d(ij+t_ldown,l,1)*ed(1)+gradq3d(ij+t_ldown,l,2)*ed(2)+gradq3d(ij+t_ldown,l,3)*ed(3) 246 475 ENDIF 247 476 qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe … … 289 518 END SUBROUTINE compute_advect_horiz 290 519 291 !==========================================================================292 PURE SUBROUTINE gradq(n0,l,n1,n2,n3,q,dq,det)293 IMPLICIT NONE294 INTEGER, INTENT(IN) :: n0,l,n1,n2,n3295 REAL(rstd), INTENT(IN) :: q(iim*jjm,llm)296 REAL(rstd), INTENT(OUT) :: dq(3), det297 298 REAL(rstd) ::detx,dety,detz299 INTEGER :: info300 INTEGER :: IPIV(3)301 302 REAL(rstd) :: A(3,3)303 REAL(rstd) :: B(3)304 305 ! TODO : replace A by A1,A2,A3306 A(1,1)=xyz_i(n1,1)-xyz_i(n0,1); A(1,2)=xyz_i(n1,2)-xyz_i(n0,2); A(1,3)=xyz_i(n1,3)-xyz_i(n0,3)307 A(2,1)=xyz_i(n2,1)-xyz_i(n0,1); A(2,2)=xyz_i(n2,2)-xyz_i(n0,2); A(2,3)=xyz_i(n2,3)-xyz_i(n0,3)308 A(3,1)=xyz_v(n3,1); A(3,2)= xyz_v(n3,2); A(3,3)=xyz_v(n3,3)309 310 dq(1) = q(n1,l)-q(n0,l)311 dq(2) = q(n2,l)-q(n0,l)312 dq(3) = 0.0313 ! CALL DGESV(3,1,A,3,IPIV,dq(:),3,info)314 CALL determinant(A(:,1),A(:,2),A(:,3),det)315 CALL determinant(dq,A(:,2),A(:,3),detx)316 CALL determinant(A(:,1),dq,A(:,3),dety)317 CALL determinant(A(:,1),A(:,2),dq,detz)318 dq(1) = detx319 dq(2) = dety320 dq(3) = detz321 END SUBROUTINE gradq322 323 !==========================================================================324 PURE SUBROUTINE determinant(a1,a2,a3,det)325 IMPLICIT NONE326 REAL(rstd), DIMENSION(3), INTENT(IN) :: a1,a2,a3327 REAL(rstd), INTENT(OUT) :: det328 REAL(rstd) :: x1,x2,x3329 x1 = a1(1) * (a2(2) * a3(3) - a2(3) * a3(2))330 x2 = a1(2) * (a2(1) * a3(3) - a2(3) * a3(1))331 x3 = a1(3) * (a2(1) * a3(2) - a2(2) * a3(1))332 det = x1 - x2 + x3333 END SUBROUTINE determinant334 520 335 521 END MODULE advect_mod -
codes/icosagcm/trunk/src/advect_tracer.f90
r174 r186 4 4 PRIVATE 5 5 6 TYPE(t_field), POINTER :: f_normal(:)7 TYPE(t_field), POINTER :: f_tangent(:)8 TYPE(t_field), POINTER :: f_gradq3d(:)9 TYPE(t_field), POINTER :: f_cc(:) ! starting point of backward-trajectory (Miura approach)10 TYPE(t_field), POINTER :: f_one_over_sqrt_leng(:)11 12 TYPE(t_message) :: req_u, req_cc, req_wfluxt, req_q, req_rhodz, req_gradq3d6 TYPE(t_field),SAVE,POINTER :: f_normal(:) 7 TYPE(t_field),SAVE,POINTER :: f_tangent(:) 8 TYPE(t_field),SAVE,POINTER :: f_gradq3d(:) 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(:) 11 12 TYPE(t_message),SAVE :: req_u, req_cc, req_wfluxt, req_q, req_rhodz, req_gradq3d 13 13 14 14 REAL(rstd), PARAMETER :: pente_max=2.0 ! for vlz 15 15 16 16 ! temporary shared variable for vlz 17 TYPE(t_field), POINTER :: f_dzqw(:) ! vertical finite difference of q18 TYPE(t_field), POINTER :: f_adzqw(:) ! abs(dzqw)19 TYPE(t_field), POINTER :: f_dzq(:) ! limited slope of q20 TYPE(t_field), POINTER :: f_wq(:) ! time-integrated flux of q17 TYPE(t_field),SAVE,POINTER :: f_dzqw(:) ! vertical finite difference of q 18 TYPE(t_field),SAVE,POINTER :: f_adzqw(:) ! abs(dzqw) 19 TYPE(t_field),SAVE,POINTER :: f_dzq(:) ! limited slope of q 20 TYPE(t_field),SAVE,POINTER :: f_wq(:) ! time-integrated flux of q 21 21 22 22 PUBLIC init_advect_tracer, advect_tracer … … 42 42 43 43 DO ind=1,ndomain 44 IF (.NOT. assigned_domain(ind)) CYCLE 44 45 CALL swap_dimensions(ind) 45 46 CALL swap_geometry(ind) … … 88 89 ENDIF 89 90 90 ! $OMP BARRIER91 !!$OMP BARRIER 91 92 92 93 CALL trace_start("advect_tracer") 93 94 94 95 CALL send_message(f_u,req_u) 96 CALL wait_message(req_u) 95 97 CALL send_message(f_wfluxt,req_wfluxt) 98 CALL wait_message(req_wfluxt) 96 99 CALL send_message(f_q,req_q) 100 CALL wait_message(req_q) 97 101 CALL send_message(f_rhodz,req_rhodz) 98 CALL wait_message(req_u)99 CALL wait_message(req_wfluxt)100 CALL wait_message(req_q)101 102 CALL wait_message(req_rhodz) 103 104 ! CALL wait_message(req_u) 105 ! CALL wait_message(req_wfluxt) 106 ! CALL wait_message(req_q) 107 ! CALL wait_message(req_rhodz) 102 108 103 109 ! 1/2 vertical transport + back-trajectories 104 110 DO ind=1,ndomain 111 IF (.NOT. assigned_domain(ind)) CYCLE 105 112 CALL swap_dimensions(ind) 106 113 CALL swap_geometry(ind) … … 126 133 127 134 CALL send_message(f_cc,req_cc) 135 CALL wait_message(req_cc) 128 136 129 137 130 138 ! horizontal transport - split in two to place transfer of gradq3d 131 ! $OMP BARRIER139 !!$OMP BARRIER 132 140 DO k = 1, nqtot 133 141 DO ind=1,ndomain 142 IF (.NOT. assigned_domain(ind)) CYCLE 134 143 CALL swap_dimensions(ind) 135 144 CALL swap_geometry(ind) … … 137 146 gradq3d = f_gradq3d(ind) 138 147 one_over_sqrt_leng=f_one_over_sqrt_leng(ind) 139 CALL compute_gradq3d(q(:,:,k),one_over_sqrt_leng,gradq3d )148 CALL compute_gradq3d(q(:,:,k),one_over_sqrt_leng,gradq3d,xyz_i,xyz_v) 140 149 END DO 141 150 142 151 CALL send_message(f_gradq3d,req_gradq3d) 143 CALL wait_message(req_cc)152 ! CALL wait_message(req_cc) 144 153 CALL wait_message(req_gradq3d) 145 154 146 155 147 156 DO ind=1,ndomain 157 IF (.NOT. assigned_domain(ind)) CYCLE 148 158 CALL swap_dimensions(ind) 149 159 CALL swap_geometry(ind) … … 158 168 159 169 ! 1/2 vertical transport 160 ! $OMP BARRIER170 !!$OMP BARRIER 161 171 162 172 DO ind=1,ndomain 173 IF (.NOT. assigned_domain(ind)) CYCLE 163 174 CALL swap_dimensions(ind) 164 175 CALL swap_geometry(ind) … … 179 190 CALL trace_end("advect_tracer") 180 191 181 ! $OMP BARRIER192 !!$OMP BARRIER 182 193 183 194 END SUBROUTINE advect_tracer … … 227 238 228 239 !--> flush dzqw, adzqw 229 ! $OMP BARRIER240 !!$OMP BARRIER 230 241 231 242 ! minmod-limited slope of q … … 260 271 261 272 !---> flush dzq 262 ! $OMP BARRIER273 !!$OMP BARRIER 263 274 264 275 ! sigw = fraction of mass that leaves level l/l+1 … … 292 303 293 304 ! --> flush wq 294 ! $OMP BARRIER305 !!$OMP BARRIER 295 306 296 307 -
codes/icosagcm/trunk/src/caldyn.f90
r162 r186 4 4 PRIVATE 5 5 CHARACTER(LEN=255),SAVE :: caldyn_type 6 !$OMP THREADPRIVATE(caldyn_type) 6 7 7 8 PUBLIC init_caldyn, caldyn, caldyn_BC -
codes/icosagcm/trunk/src/caldyn_adv.f90
r146 r186 27 27 28 28 DO ind=1,ndomain 29 IF (.NOT. assigned_domain(ind)) CYCLE 29 30 CALL swap_dimensions(ind) 30 31 CALL swap_geometry(ind) … … 78 79 79 80 DO ind=1,ndomain 81 IF (.NOT. assigned_domain(ind)) CYCLE 80 82 CALL swap_dimensions(ind) 81 83 CALL swap_geometry(ind) … … 88 90 du=f_du(ind) 89 91 90 !$OMP PARALLEL DEFAULT(SHARED)92 ! !$OMP PARALLEL DEFAULT(SHARED) 91 93 CALL compute_caldyn(ps,u,hflux, wflux, dps, dtheta_rhodz, du) 92 !$OMP END PARALLEL94 ! !$OMP END PARALLEL 93 95 ENDDO 94 96 -
codes/icosagcm/trunk/src/caldyn_gcm.f90
r183 r186 6 6 INTEGER, PARAMETER :: energy=1, enstrophy=2 7 7 TYPE(t_field),POINTER :: f_out_u(:), f_qu(:), f_qv(:) 8 REAL(rstd),POINTER :: out_u(:,:), p(:,:), qu(:,:) 8 REAL(rstd),SAVE,POINTER :: out_u(:,:), p(:,:), qu(:,:) 9 !$OMP THREADPRIVATE(out_u, p, qu) 9 10 10 11 TYPE(t_field),POINTER :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:) … … 19 20 20 21 INTEGER :: caldyn_conserv 21 22 !$OMP THREADPRIVATE(caldyn_conserv) 23 22 24 TYPE(t_message) :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu 23 25 … … 101 103 IF (omp_first) THEN 102 104 DO ind=1,ndomain 105 IF (.NOT. assigned_domain(ind)) CYCLE 103 106 CALL swap_dimensions(ind) 104 107 CALL swap_geometry(ind) … … 124 127 ENDIF 125 128 126 !$OMP BARRIER129 ! !$OMP BARRIER 127 130 END SUBROUTINE caldyn_BC 128 131 … … 175 178 IF(caldyn_eta==eta_mass) THEN 176 179 CALL init_message(f_ps,req_i1,req_ps) 177 ELSE180 ! ELSE 178 181 CALL init_message(f_mass,req_i1,req_mass) 179 182 END IF … … 181 184 CALL init_message(f_u,req_e1_vect,req_u) 182 185 CALL init_message(f_qu,req_e1_scal,req_qu) 186 ! IF(caldyn_eta==eta_mass) THEN 187 ! CALL send_message(f_ps,req_ps) 188 ! CALL wait_message(req_ps) 189 ! ELSE 190 ! CALL send_message(f_mass,req_mass) 191 ! CALL wait_message(req_mass) 192 ! END IF 193 ENDIF 194 195 CALL trace_start("caldyn") 196 183 197 IF(caldyn_eta==eta_mass) THEN 184 198 CALL send_message(f_ps,req_ps) 199 CALL wait_message(req_ps) 185 200 ELSE 186 201 CALL send_message(f_mass,req_mass) 202 CALL wait_message(req_mass) 187 203 END IF 188 ENDIF189 190 CALL trace_start("caldyn")191 204 192 205 CALL send_message(f_u,req_u) 206 CALL wait_message(req_u) 193 207 CALL send_message(f_theta_rhodz,req_theta_rhodz) 208 CALL wait_message(req_theta_rhodz) 209 210 ! CALL wait_message(req_u) 211 ! CALL wait_message(req_theta_rhodz) 194 212 195 213 SELECT CASE(caldyn_conserv) 196 214 CASE(energy) ! energy-conserving 197 215 DO ind=1,ndomain 216 IF (.NOT. assigned_domain(ind)) CYCLE 198 217 CALL swap_dimensions(ind) 199 218 CALL swap_geometry(ind) … … 209 228 210 229 CALL send_message(f_qu,req_qu) 230 CALL wait_message(req_qu) 211 231 212 232 DO ind=1,ndomain 233 IF (.NOT. assigned_domain(ind)) CYCLE 213 234 CALL swap_dimensions(ind) 214 235 CALL swap_geometry(ind) … … 237 258 CASE(enstrophy) ! enstrophy-conserving 238 259 DO ind=1,ndomain 260 IF (.NOT. assigned_domain(ind)) CYCLE 239 261 CALL swap_dimensions(ind) 240 262 CALL swap_geometry(ind) … … 267 289 END SELECT 268 290 269 ! $OMP BARRIER291 !!$OMP BARRIER 270 292 IF (write_out) THEN 271 293 … … 289 311 ! CALL check_mass_conservation(f_ps,f_dps) 290 312 CALL trace_end("caldyn") 291 ! $OMP BARRIER313 !!$OMP BARRIER 292 314 293 315 END SUBROUTINE caldyn … … 338 360 339 361 IF(caldyn_eta==eta_mass) THEN 340 CALL wait_message(req_ps)362 ! CALL wait_message(req_ps) 341 363 ELSE 342 CALL wait_message(req_mass)364 ! CALL wait_message(req_mass) 343 365 END IF 344 CALL wait_message(req_theta_rhodz)366 ! CALL wait_message(req_theta_rhodz) 345 367 346 368 IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta 347 369 DO l = ll_begin,ll_end 348 CALL test_message(req_u)370 ! CALL test_message(req_u) 349 371 !DIR$ SIMD 350 372 DO ij=ij_begin_ext,ij_end_ext … … 356 378 ELSE ! Compute only theta 357 379 DO l = ll_begin,ll_end 358 CALL test_message(req_u)380 ! CALL test_message(req_u) 359 381 !DIR$ SIMD 360 382 DO ij=ij_begin_ext,ij_end_ext … … 364 386 END IF 365 387 366 CALL wait_message(req_u)388 ! CALL wait_message(req_u) 367 389 368 390 !!! Compute shallow-water potential vorticity … … 427 449 !!! Compute exner function and geopotential 428 450 DO l = 1,llm 429 451 ! !$OMP DO SCHEDULE(STATIC) 430 452 !DIR$ SIMD 431 453 DO ij=ij_begin_ext,ij_end_ext … … 448 470 ! specific volume 1 = dphi/g/rhodz 449 471 DO l = 1,llm 450 !$OMP DO SCHEDULE(STATIC)472 ! !$OMP DO SCHEDULE(STATIC) 451 473 !DIR$ SIMD 452 474 DO ij=ij_begin_ext,ij_end_ext … … 462 484 ! other layers 463 485 DO l = llm-1, 1, -1 464 !$OMP DO SCHEDULE(STATIC) 486 487 ! !$OMP DO SCHEDULE(STATIC) 465 488 !DIR$ SIMD 466 489 DO ij=ij_begin_ext,ij_end_ext … … 475 498 ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 476 499 DO l = 1,llm 477 !$OMP DO SCHEDULE(STATIC) 500 501 ! !$OMP DO SCHEDULE(STATIC) 478 502 !DIR$ SIMD 479 503 DO ij=ij_begin_ext,ij_end_ext … … 521 545 CALL trace_start("compute_caldyn_horiz") 522 546 523 CALL wait_message(req_theta_rhodz)547 ! CALL wait_message(req_theta_rhodz) 524 548 525 549 DO l = ll_begin, ll_end 526 550 !!! Compute mass and theta fluxes 527 IF (caldyn_conserv==energy) CALL test_message(req_qu)551 ! IF (caldyn_conserv==energy) CALL test_message(req_qu) 528 552 !DIR$ SIMD 529 553 DO ij=ij_begin_ext,ij_end_ext … … 565 589 CASE(energy) ! energy-conserving TRiSK 566 590 567 CALL wait_message(req_qu)591 ! CALL wait_message(req_qu) 568 592 569 593 DO l=ll_begin,ll_end … … 671 695 ! other layers 672 696 DO l = llm-1, 1, -1 673 !$OMP DO SCHEDULE(STATIC)697 ! !$OMP DO SCHEDULE(STATIC) 674 698 !DIR$ SIMD 675 699 DO ij=ij_begin_ext,ij_end_ext … … 686 710 !DIR$ SIMD 687 711 DO ij=ij_begin,ij_end 688 berni(ij,l) = pk(ij,l) & 689 + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 + & 712 713 berni(ij,l) = pk(ij,l) + & 714 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 + & 690 715 le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 + & 691 716 le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 + & … … 756 781 REAL(rstd),INTENT(INOUT) :: convm(iim*jjm,llm) ! mass flux convergence 757 782 758 REAL(rstd),INTENT( OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)759 REAL(rstd),INTENT( OUT) :: wwuu(iim*3*jjm,llm+1)783 REAL(rstd),INTENT(INOUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s) 784 REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1) 760 785 REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm) 761 786 REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm) … … 766 791 REAL(rstd) :: p_ik, exner_ik 767 792 768 793 ! REAL(rstd) :: wwuu(iim*3*jjm,llm+1) ! tmp var, don't know why but gain 30% on the whole code in opemp 794 ! need to be understood 795 796 ! wwuu=wwuu_out 769 797 CALL trace_start("compute_caldyn_vert") 770 798 771 ! $OMP BARRIER799 !!$OMP BARRIER 772 800 !!! cumulate mass flux convergence from top to bottom 773 801 DO l = llm-1, 1, -1 774 IF (caldyn_conserv==energy) CALL test_message(req_qu) 775 !$OMP DO SCHEDULE(STATIC) 802 ! IF (caldyn_conserv==energy) CALL test_message(req_qu) 803 804 !!$OMP DO SCHEDULE(STATIC) 776 805 !DIR$ SIMD 777 806 DO ij=ij_begin,ij_end … … 794 823 !!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC) 795 824 DO l=ll_beginp1,ll_end 796 IF (caldyn_conserv==energy) CALL test_message(req_qu)825 ! IF (caldyn_conserv==energy) CALL test_message(req_qu) 797 826 !DIR$ SIMD 798 827 DO ij=ij_begin,ij_end … … 828 857 829 858 !--> flush wwuu 830 ! $OMP BARRIER859 ! !$OMP BARRIER 831 860 832 861 ! Add vertical transport to du … … 839 868 ENDDO 840 869 ENDDO 870 871 ! DO l=ll_beginp1,ll_end 872 !!DIR$ SIMD 873 ! DO ij=ij_begin,ij_end 874 ! wwuu_out(ij+u_right,l) = wwuu(ij+u_right,l) 875 ! wwuu_out(ij+u_lup,l) = wwuu(ij+u_lup,l) 876 ! wwuu_out(ij+u_ldown,l) = wwuu(ij+u_ldown,l) 877 ! ENDDO 878 ! ENDDO 841 879 842 880 CALL trace_end("compute_caldyn_vert") … … 973 1011 974 1012 DO ind=1,ndomain 1013 IF (.NOT. assigned_domain(ind)) CYCLE 975 1014 CALL swap_dimensions(ind) 976 1015 CALL swap_geometry(ind) … … 1002 1041 1003 1042 DO ind=1,ndomain 1043 IF (.NOT. assigned_domain(ind)) CYCLE 1004 1044 CALL swap_dimensions(ind) 1005 1045 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/check_conserve.f90
r172 r186 10 10 PUBLIC init_check_conserve, check_conserve 11 11 REAL(rstd),SAVE:: mtot0,ztot0,etot0,ang0,stot0,rmsv0 12 !$OMP THREADPRIVATE(mtot0,ztot0,etot0,ang0,stot0,rmsv0) 12 13 REAL(rstd),SAVE:: etot,ang,stot,rmsv 14 !$OMP THREADPRIVATE(etot,ang,stot,rmsv) 13 15 REAL(rstd),SAVE:: ztot 14 16 !$OMP THREADPRIVATE(ztot) 15 17 16 18 CONTAINS … … 53 55 54 56 DO ind=1,ndomain 57 IF (.NOT. assigned_domain(ind)) CYCLE 55 58 CALL swap_dimensions(ind) 56 59 CALL swap_geometry(ind) … … 67 70 68 71 IF (is_mpi_root) THEN 69 72 !$OMP MASTER 70 73 IF ( it == 0 ) Then 71 74 ztot0 = ztot … … 87 90 WRITE(*,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang 88 91 89 4000 FORMAT(10x,'masse', 4x,'rmsdpdt',7x,'energie',2x,'enstrophie' &90 , 2x,'entropie',3x,'rmsv',4x,'mt.ang',/,'GLOB ' &92 4000 FORMAT(10x,'masse',5x,'rmsdpdt',5x,'energie',5x,'enstrophie' & 93 ,5x,'entropie',5x,'rmsv',5x,'mt.ang',/,'GLOB ' & 91 94 ,e10.3,e13.6,5e10.3/) 92 close(134) 95 close(134) 96 !$OMP END MASTER 93 97 END IF 94 98 END SUBROUTINE check_conserve … … 99 103 USE mpi_mod 100 104 USE mpipara 105 USE transfert_omp_mod 101 106 USE icosa 102 107 IMPLICIT NONE … … 108 113 INTEGER :: ind,i,j,ij 109 114 REAL :: mloc, rmsloc 115 REAL :: mloc_mpi, rmsloc_mpi 110 116 111 117 mloc=0.0; rmsloc=0.0 112 118 DO ind=1,ndomain 119 IF (.NOT. assigned_domain(ind)) CYCLE 113 120 CALL swap_dimensions(ind) 114 121 CALL swap_geometry(ind) … … 127 134 128 135 IF (using_mpi) THEN 129 CALL MPI_REDUCE(mloc, mtot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 130 CALL MPI_REDUCE(rmsloc, rmsdpdt, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 136 CALL reduce_sum_omp(mloc, mloc_mpi) 137 CALL reduce_sum_omp(rmsloc, rmsloc_mpi) 138 !$OMP MASTER 139 CALL MPI_REDUCE(mloc_mpi, mtot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 140 CALL MPI_REDUCE(rmsloc_mpi, rmsdpdt, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 141 !$OMP END MASTER 142 ELSE 143 131 144 ENDIF 132 145 … … 153 166 154 167 DO ind=1,ndomain 168 IF (.NOT. assigned_domain(ind)) CYCLE 155 169 CALL swap_dimensions(ind) 156 170 CALL swap_geometry(ind) … … 260 274 261 275 DO ind=1,ndomain 276 IF (.NOT. assigned_domain(ind)) CYCLE 262 277 CALL swap_dimensions(ind) 263 278 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/dcmip_initial_conditions_test_1_2_3_v5.f90
r54 r186 260 260 r = ACOS (sin_tmp + cos_tmp*cos(lon-lambda0)) 261 261 r2 = ACOS (sin_tmp2 + cos_tmp2*cos(lon-lambda1)) 262 d1 = min( 1. d0, (r/RR)**2 + ((height-z0)/ZZ)**2 )263 d2 = min( 1. d0, (r2/RR)**2 + ((height-z0)/ZZ)**2 )262 d1 = min( 1., (r/RR)**2 + ((height-z0)/ZZ)**2 ) 263 d2 = min( 1., (r2/RR)**2 + ((height-z0)/ZZ)**2 ) 264 264 265 265 q1 = 0.5d0 * (1.d0 + cos(pi*d1)) + 0.5d0 * (1.d0 + cos(pi*d2)) -
codes/icosagcm/trunk/src/dimensions.f90
r174 r186 1 1 MODULE dimensions 2 2 3 INTEGER :: iim 4 INTEGER :: jjm 5 INTEGER :: ii_begin 6 INTEGER :: jj_begin 7 INTEGER :: ii_end 8 INTEGER :: jj_end 9 INTEGER :: ii_nb 10 INTEGER :: jj_nb 11 INTEGER :: ij_begin 12 INTEGER :: ij_end 13 INTEGER :: ij_begin_ext 14 INTEGER :: ij_end_ext 3 INTEGER,SAVE :: iim 4 !$OMP THREADPRIVATE(iim) 5 INTEGER,SAVE :: jjm 6 !$OMP THREADPRIVATE(jjm) 7 INTEGER,SAVE :: ii_begin 8 !$OMP THREADPRIVATE(ii_begin) 9 INTEGER,SAVE :: jj_begin 10 !$OMP THREADPRIVATE(jj_begin) 11 INTEGER,SAVE :: ii_end 12 !$OMP THREADPRIVATE(ii_end) 13 INTEGER,SAVE :: jj_end 14 !$OMP THREADPRIVATE(jj_end) 15 INTEGER,SAVE :: ii_nb 16 !$OMP THREADPRIVATE(ii_nb) 17 INTEGER,SAVE :: jj_nb 18 !$OMP THREADPRIVATE(jj_nb) 19 INTEGER,SAVE :: ij_begin 20 !$OMP THREADPRIVATE(ij_begin) 21 INTEGER,SAVE :: ij_end 22 !$OMP THREADPRIVATE(ij_end) 23 INTEGER,SAVE :: ij_begin_ext 24 !$OMP THREADPRIVATE(ij_begin_ext) 25 INTEGER,SAVE :: ij_end_ext 26 !$OMP THREADPRIVATE(ij_end_ext) 15 27 16 INTEGER :: t_right 17 INTEGER :: t_rup 18 INTEGER :: t_lup 19 INTEGER :: t_left 20 INTEGER :: t_ldown 21 INTEGER :: t_rdown 28 INTEGER,SAVE :: t_right 29 !$OMP THREADPRIVATE(t_right) 30 INTEGER,SAVE :: t_rup 31 !$OMP THREADPRIVATE(t_rup) 32 INTEGER,SAVE :: t_lup 33 !$OMP THREADPRIVATE(t_lup) 34 INTEGER,SAVE :: t_left 35 !$OMP THREADPRIVATE(t_left) 36 INTEGER,SAVE :: t_ldown 37 !$OMP THREADPRIVATE(t_ldown) 38 INTEGER,SAVE :: t_rdown 39 !$OMP THREADPRIVATE(t_rdown) 22 40 23 INTEGER :: u_right 24 INTEGER :: u_rup 25 INTEGER :: u_lup 26 INTEGER :: u_left 27 INTEGER :: u_ldown 28 INTEGER :: u_rdown 41 INTEGER,SAVE :: u_right 42 !$OMP THREADPRIVATE(u_right) 43 INTEGER,SAVE :: u_rup 44 !$OMP THREADPRIVATE(u_rup) 45 INTEGER,SAVE :: u_lup 46 !$OMP THREADPRIVATE(u_lup) 47 INTEGER,SAVE :: u_left 48 !$OMP THREADPRIVATE(u_left) 49 INTEGER,SAVE :: u_ldown 50 !$OMP THREADPRIVATE(u_ldown) 51 INTEGER,SAVE :: u_rdown 52 !$OMP THREADPRIVATE(u_rdown) 29 53 30 INTEGER :: z_rup 31 INTEGER :: z_up 32 INTEGER :: z_lup 33 INTEGER :: z_ldown 34 INTEGER :: z_down 35 INTEGER :: z_rdown 54 INTEGER,SAVE :: z_rup 55 !$OMP THREADPRIVATE(z_rup) 56 INTEGER,SAVE :: z_up 57 !$OMP THREADPRIVATE(z_up) 58 INTEGER,SAVE :: z_lup 59 !$OMP THREADPRIVATE(z_lup) 60 INTEGER,SAVE :: z_ldown 61 !$OMP THREADPRIVATE(z_ldown) 62 INTEGER,SAVE :: z_down 63 !$OMP THREADPRIVATE(z_down) 64 INTEGER,SAVE :: z_rdown 65 !$OMP THREADPRIVATE(z_rdown) 36 66 37 INTEGER :: t_pos(6) 38 INTEGER :: u_pos(6) 39 INTEGER :: z_pos(6) 67 INTEGER,SAVE :: t_pos(6) 68 !$OMP THREADPRIVATE(t_pos) 69 INTEGER,SAVE :: u_pos(6) 70 !$OMP THREADPRIVATE(u_pos) 71 INTEGER,SAVE :: z_pos(6) 72 !$OMP THREADPRIVATE(z_pos) 40 73 41 74 CONTAINS … … 48 81 TYPE(t_domain),POINTER :: d 49 82 50 51 !$OMP MASTER52 83 d=>domain(ind) 53 84 … … 91 122 z_pos(:)=d%z_pos(:) 92 123 93 !$OMP END MASTER94 !$OMP BARRIER95 124 END SUBROUTINE swap_dimensions 96 125 -
codes/icosagcm/trunk/src/dissip_gcm.f90
r175 r186 10 10 TYPE(t_field),POINTER,SAVE :: f_dtheta_diss(:) 11 11 TYPE(t_field),POINTER,SAVE :: f_dtheta_rhodz_diss(:) 12 TYPE(t_message) :: req_due, req_dtheta12 TYPE(t_message),SAVE :: req_due, req_dtheta 13 13 14 14 INTEGER,SAVE :: nitergdiv=1 15 !$OMP THREADPRIVATE(nitergdiv) 15 16 INTEGER,SAVE :: nitergrot=1 17 !$OMP THREADPRIVATE(nitergrot) 16 18 INTEGER,SAVE :: niterdivgrad=1 19 !$OMP THREADPRIVATE(niterdivgrad) 17 20 18 21 REAL,ALLOCATABLE,SAVE :: tau_graddiv(:) 22 !$OMP THREADPRIVATE(tau_graddiv) 19 23 REAL,ALLOCATABLE,SAVE :: tau_gradrot(:) 24 !$OMP THREADPRIVATE(tau_gradrot) 20 25 REAL,ALLOCATABLE,SAVE :: tau_divgrad(:) 26 !$OMP THREADPRIVATE(tau_divgrad) 21 27 22 28 REAL,SAVE :: cgraddiv 29 !$OMP THREADPRIVATE(cgraddiv) 23 30 REAL,SAVE :: cgradrot 31 !$OMP THREADPRIVATE(cgradrot) 24 32 REAL,SAVE :: cdivgrad 33 !$OMP THREADPRIVATE(cdivgrad) 25 34 26 35 INTEGER, SAVE :: rayleigh_friction_type, rayleigh_shear 27 REAL, save :: rayleigh_tau 36 !$OMP THREADPRIVATE(rayleigh_friction_type) 37 REAL, SAVE :: rayleigh_tau 38 !$OMP THREADPRIVATE(rayleigh_shear) 28 39 29 40 REAL,SAVE :: dtdissip 41 !$OMP THREADPRIVATE(dtdissip) 30 42 31 43 PUBLIC init_dissip, dissip … … 58 70 USE transfert_mod 59 71 USE time_mod 60 72 USE transfert_omp_mod 61 73 IMPLICIT NONE 62 74 63 TYPE(t_field),POINTER :: f_u(:)64 TYPE(t_field),POINTER :: f_du(:)65 REAL(rstd),POINTER :: u(:)66 REAL(rstd),POINTER :: du(:)67 TYPE(t_field),POINTER :: f_theta(:)68 TYPE(t_field),POINTER :: f_dtheta(:)75 TYPE(t_field),POINTER,SAVE :: f_u(:) 76 TYPE(t_field),POINTER,SAVE :: f_du(:) 77 REAL(rstd),POINTER :: u(:) 78 REAL(rstd),POINTER :: du(:) 79 TYPE(t_field),POINTER,SAVE :: f_theta(:) 80 TYPE(t_field),POINTER ,SAVE :: f_dtheta(:) 69 81 REAL(rstd),POINTER :: theta(:) 70 82 REAL(rstd),POINTER :: dtheta(:) … … 128 140 129 141 tau_gradrot(:)=5000 130 CALL getin("tau_gradrot",tau _gradrot)142 CALL getin("tau_gradrot",tau) 131 143 tau_gradrot(:)=tau/scale_factor 132 144 … … 146 158 147 159 DO ind=1,ndomain 160 IF (.NOT. assigned_domain(ind)) CYCLE 148 161 CALL swap_dimensions(ind) 149 162 CALL swap_geometry(ind) … … 171 184 CALL transfert_request(f_u,req_e1_vect) 172 185 DO ind=1,ndomain 186 IF (.NOT. assigned_domain(ind)) CYCLE 173 187 CALL swap_dimensions(ind) 174 188 CALL swap_geometry(ind) … … 183 197 184 198 DO ind=1,ndomain 199 IF (.NOT. assigned_domain(ind)) CYCLE 185 200 CALL swap_dimensions(ind) 186 201 CALL swap_geometry(ind) … … 199 214 200 215 IF (using_mpi) THEN 201 CALL MPI_ALLREDUCE(dumax,dumax1,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 216 CALL reduce_sum_omp(dumax,dumax1) 217 !$OMP MASTER 218 CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 219 !$OMP END MASTER 220 CALL bcast_omp(dumax) 221 ELSE 222 CALL allreduce_sum_omp(dumax,dumax1) 202 223 dumax=dumax1 203 ENDIF 224 ENDIF 204 225 205 226 DO ind=1,ndomain 227 IF (.NOT. assigned_domain(ind)) CYCLE 206 228 CALL swap_dimensions(ind) 207 229 CALL swap_geometry(ind) … … 223 245 224 246 DO ind=1,ndomain 247 IF (.NOT. assigned_domain(ind)) CYCLE 225 248 CALL swap_dimensions(ind) 226 249 CALL swap_geometry(ind) … … 247 270 CALL transfert_request(f_u,req_e1_vect) 248 271 DO ind=1,ndomain 272 IF (.NOT. assigned_domain(ind)) CYCLE 249 273 CALL swap_dimensions(ind) 250 274 CALL swap_geometry(ind) … … 259 283 260 284 DO ind=1,ndomain 285 IF (.NOT. assigned_domain(ind)) CYCLE 261 286 CALL swap_dimensions(ind) 262 287 CALL swap_geometry(ind) … … 274 299 ENDDO 275 300 276 IF (using_mpi) THEN 277 CALL MPI_ALLREDUCE(dumax,dumax1,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 301 IF (using_mpi) THEN 302 CALL reduce_sum_omp(dumax,dumax1) 303 !$OMP MASTER 304 CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 305 !$OMP END MASTER 306 CALL bcast_omp(dumax) 307 ELSE 308 CALL allreduce_sum_omp(dumax,dumax1) 278 309 dumax=dumax1 279 ENDIF 310 ENDIF 311 280 312 281 313 DO ind=1,ndomain 314 IF (.NOT. assigned_domain(ind)) CYCLE 282 315 CALL swap_dimensions(ind) 283 316 CALL swap_geometry(ind) … … 298 331 299 332 DO ind=1,ndomain 333 IF (.NOT. assigned_domain(ind)) CYCLE 300 334 CALL swap_dimensions(ind) 301 335 CALL swap_geometry(ind) … … 317 351 CALL transfert_request(f_theta,req_i1) 318 352 DO ind=1,ndomain 353 IF (.NOT. assigned_domain(ind)) CYCLE 319 354 CALL swap_dimensions(ind) 320 355 CALL swap_geometry(ind) … … 329 364 330 365 DO ind=1,ndomain 366 IF (.NOT. assigned_domain(ind)) CYCLE 331 367 CALL swap_dimensions(ind) 332 368 CALL swap_geometry(ind) … … 341 377 ENDDO 342 378 ENDDO 343 344 IF (using_mpi) THEN 345 CALL MPI_ALLREDUCE(dthetamax,dthetamax1,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 346 dthetamax=dthetamax1 347 ENDIF 379 380 IF (using_mpi) THEN 381 CALL reduce_sum_omp(dthetamax ,dthetamax1) 382 !$OMP MASTER 383 CALL MPI_ALLREDUCE(dthetamax1,dthetamax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 384 !$OMP END MASTER 385 CALL bcast_omp(dthetamax) 386 ELSE 387 CALL allreduce_sum_omp(dthetamax,dthetamax1) 388 dumax=dumax1 389 ENDIF 348 390 349 391 IF (is_mpi_root) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax 350 392 351 393 DO ind=1,ndomain 394 IF (.NOT. assigned_domain(ind)) CYCLE 352 395 CALL swap_dimensions(ind) 353 396 CALL swap_geometry(ind) … … 434 477 435 478 DO ind=1,ndomain 479 IF (.NOT. assigned_domain(ind)) CYCLE 436 480 CALL swap_dimensions(ind) 437 481 CALL swap_geometry(ind) … … 524 568 525 569 DO ind=1,ndomain 570 IF (.NOT. assigned_domain(ind)) CYCLE 526 571 CALL swap_dimensions(ind) 527 572 CALL swap_geometry(ind) … … 544 589 545 590 DO ind=1,ndomain 591 IF (.NOT. assigned_domain(ind)) CYCLE 546 592 CALL swap_dimensions(ind) 547 593 CALL swap_geometry(ind) … … 571 617 572 618 DO ind=1,ndomain 619 IF (.NOT. assigned_domain(ind)) CYCLE 573 620 CALL swap_dimensions(ind) 574 621 CALL swap_geometry(ind) … … 591 638 592 639 DO ind=1,ndomain 640 IF (.NOT. assigned_domain(ind)) CYCLE 593 641 CALL swap_dimensions(ind) 594 642 CALL swap_geometry(ind) … … 618 666 619 667 DO ind=1,ndomain 668 IF (.NOT. assigned_domain(ind)) CYCLE 620 669 CALL swap_dimensions(ind) 621 670 CALL swap_geometry(ind) … … 630 679 631 680 DO ind=1,ndomain 681 IF (.NOT. assigned_domain(ind)) CYCLE 632 682 CALL swap_dimensions(ind) 633 683 CALL swap_geometry(ind) … … 661 711 662 712 DO ind=1,ndomain 713 IF (.NOT. assigned_domain(ind)) CYCLE 663 714 CALL swap_dimensions(ind) 664 715 CALL swap_geometry(ind) … … 679 730 CALL wait_message(req_dtheta) 680 731 DO ind=1,ndomain 732 IF (.NOT. assigned_domain(ind)) CYCLE 681 733 CALL swap_dimensions(ind) 682 734 CALL swap_geometry(ind) … … 688 740 689 741 DO ind=1,ndomain 742 IF (.NOT. assigned_domain(ind)) CYCLE 690 743 CALL swap_dimensions(ind) 691 744 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/disvert.f90
r177 r186 2 2 USE icosa 3 3 REAL(rstd), SAVE, POINTER :: ap(:) 4 !$OMP THREADPRIVATE(ap) 4 5 REAL(rstd), SAVE, POINTER :: bp(:) 6 !$OMP THREADPRIVATE(bp) 5 7 REAL(rstd), SAVE, POINTER :: presnivs(:) 8 !$OMP THREADPRIVATE(presnivs) 6 9 REAL(rstd), SAVE, POINTER :: mass_al(:), mass_bl(:), mass_ak(:), mass_bk(:), mass_dak(:), mass_dbk(:) 10 !$OMP THREADPRIVATE(mass_al, mass_bl, mass_ak, mass_bk, mass_dak, mass_dbk) 7 11 8 12 REAL(rstd) :: ptop ! pressure at top of atmosphere l=llm+1 13 !$OMP THREADPRIVATE(ptop) 9 14 ! vertical coordinate : Lagrangian or mass-based 10 15 INTEGER, SAVE :: caldyn_eta 16 !$OMP THREADPRIVATE(caldyn_eta) 11 17 INTEGER, PARAMETER :: eta_mass=1, eta_lag=2 12 18 LOGICAL,SAVE :: ap_bp_present 19 !$OMP THREADPRIVATE(ap_bp_present) 13 20 14 21 CONTAINS … … 42 49 ap_bp_present=.TRUE. 43 50 CALL getin("disvert",def) 51 PRINT*,"def --> ",def 44 52 SELECT CASE (TRIM(def)) 45 53 CASE('none') … … 83 91 bp=>bp_ncarl30 84 92 presnivs=>presnivs_ncarl30 85 93 PRINT*,"ap --> ",ap 94 PRINT*,"ap_ncarl30 --> ",ap_ncarl30 86 95 CASE default 87 96 IF (is_mpi_root) PRINT*,'Bad selector for variable disvert : <', TRIM(def),"> options are <std>, <ncar>, <ncarl30>" … … 169 178 INTEGER :: ncid,levid,ilevid,hyaiid,hybiid,hyamid,hybmid,P0id 170 179 INTEGER :: l 171 172 180 181 !$OMP BARRIER 182 !$OMP MASTER 183 IF(ap_bp_present) THEN 173 184 status = NF90_CREATE('apbp.nc', NF90_CLOBBER, ncid) 174 185 status = NF90_DEF_DIM(ncid,'lev',llm,lev) … … 204 215 205 216 status = NF90_ENDDEF(ncid) 206 207 IF(ap_bp_present) THEN 217 208 218 status=NF90_PUT_VAR(ncid,ilevid, ap(:)+bp(:)*Preff) 209 219 … … 228 238 229 239 status=NF90_PUT_VAR(ncid,P0id, Preff) 230 231 END IF232 240 233 241 status=NF90_CLOSE(ncid) 234 242 END IF 243 !$OMP END MASTER 244 !$OMP BARRIER 235 245 END SUBROUTINE write_apbp 236 246 -
codes/icosagcm/trunk/src/disvert_dcmip200.f90
r131 r186 3 3 4 4 REAL(rstd), SAVE, ALLOCATABLE,TARGET :: ap(:) 5 !$OMP THREADPRIVATE(ap) 5 6 REAL(rstd), SAVE, ALLOCATABLE,TARGET :: bp(:) 7 !$OMP THREADPRIVATE(bp) 6 8 REAL(rstd), SAVE, ALLOCATABLE,TARGET :: presnivs(:) 9 !$OMP THREADPRIVATE(presnivs) 7 10 8 11 CONTAINS -
codes/icosagcm/trunk/src/disvert_dcmip3.f90
r131 r186 3 3 4 4 REAL(rstd), SAVE, ALLOCATABLE,TARGET :: ap(:) 5 !$OMP THREADPRIVATE(ap) 5 6 REAL(rstd), SAVE, ALLOCATABLE,TARGET :: bp(:) 7 !$OMP THREADPRIVATE(bp) 6 8 REAL(rstd), SAVE, ALLOCATABLE,TARGET :: presnivs(:) 9 !$OMP THREADPRIVATE(presnivs) 7 10 8 11 CONTAINS -
codes/icosagcm/trunk/src/disvert_ncar.f90
r131 r186 3 3 4 4 REAL(rstd), SAVE, ALLOCATABLE,TARGET :: ap(:) 5 !$OMP THREADPRIVATE(ap) 5 6 REAL(rstd), SAVE, ALLOCATABLE,TARGET :: bp(:) 7 !$OMP THREADPRIVATE(bp) 6 8 REAL(rstd), SAVE, ALLOCATABLE,TARGET :: presnivs(:) 9 !$OMP THREADPRIVATE(presnivs) 7 10 8 11 CONTAINS -
codes/icosagcm/trunk/src/disvert_ncarl30.f90
r131 r186 1 1 MODULE disvert_ncarl30_mod 2 2 USE icosa 3 3 4 REAL(rstd), SAVE, ALLOCATABLE,TARGET :: ap(:) 5 !$OMP THREADPRIVATE(ap) 4 6 REAL(rstd), SAVE, ALLOCATABLE,TARGET :: bp(:) 7 !$OMP THREADPRIVATE(bp) 5 8 REAL(rstd), SAVE, ALLOCATABLE,TARGET :: presnivs(:) 9 !$OMP THREADPRIVATE(presnivs) 6 10 7 11 CONTAINS -
codes/icosagcm/trunk/src/disvert_std.f90
r131 r186 2 2 USE icosa 3 3 REAL(rstd), SAVE, ALLOCATABLE,TARGET :: ap(:) 4 !$OMP THREADPRIVATE(ap) 4 5 REAL(rstd), SAVE, ALLOCATABLE,TARGET :: bp(:) 6 !$OMP THREADPRIVATE(bp) 5 7 REAL(rstd), SAVE, ALLOCATABLE,TARGET :: presnivs(:) 8 !$OMP THREADPRIVATE(presnivs) 6 9 7 10 CONTAINS -
codes/icosagcm/trunk/src/domain.f90
r174 r186 58 58 END TYPE t_domain 59 59 60 INTEGER,SAVE :: ndomain 61 INTEGER,SAVE :: ndomain_glo 62 60 63 TYPE(t_domain),SAVE,ALLOCATABLE,TARGET :: domain(:) 61 INTEGER :: ndomain62 64 TYPE(t_domain),SAVE,ALLOCATABLE,TARGET :: domain_glo(:) 63 INTEGER :: ndomain_glo 64 65 INTEGER,ALLOCATABLE,SAVE :: domglo_rank(:) 66 INTEGER,ALLOCATABLE,SAVE :: domglo_loc_ind(:) 67 INTEGER,ALLOCATABLE,SAVE :: domloc_glo_ind(:) 68 65 66 INTEGER,ALLOCATABLE,SAVE :: domglo_rank(:) ! size : ndomain_glo : mpi rank assigned to a domain 67 INTEGER,ALLOCATABLE,SAVE :: domglo_loc_ind(:) ! size : ndomain_glo : corresponding local indice for a global domain indice 68 INTEGER,ALLOCATABLE,SAVE :: domloc_glo_ind(:) ! size : ndomain : corresponding global indice for a local domain indice 69 70 LOGICAL, ALLOCATABLE, SAVE :: assigned_domain(:) ! size : ndomain : is an omp task is assigned to this domain 71 !$OMP THREADPRIVATE(assigned_domain) 72 69 73 CONTAINS 70 74 … … 423 427 ENDIF 424 428 ENDDO 429 430 !$OMP PARALLEL 431 CALL assign_domain_omp 432 !$OMP END PARALLEL 433 425 434 426 435 END SUBROUTINE assign_domain 427 436 437 SUBROUTINE assign_domain_omp 438 USE omp_para 439 USE mpipara 440 IMPLICIT NONE 441 INTEGER :: nb_domain 442 INTEGER :: rank, ind, i 443 444 ALLOCATE(assigned_domain(ndomain)) 445 446 ind=0 447 DO rank=0,omp_size-1 448 nb_domain=ndomain/omp_size 449 IF ( rank < MOD(ndomain,omp_size) ) nb_domain=nb_domain+1 450 451 DO i=1,nb_domain 452 ind=ind+1 453 IF (rank==omp_rank) THEN 454 assigned_domain(ind)=.TRUE. 455 PRINT *,"Rank ",mpi_rank," task ",rank," local domain : ",ind," global domain ",domloc_glo_ind(ind) 456 ELSE 457 assigned_domain(ind)=.FALSE. 458 ENDIF 459 ENDDO 460 461 ENDDO 462 463 END SUBROUTINE assign_domain_omp 464 465 428 466 429 467 SUBROUTINE compute_domain -
codes/icosagcm/trunk/src/earth_const.f90
r165 r186 18 18 19 19 SUBROUTINE init_earth_const 20 USE ioipsl20 USE getin_mod 21 21 IMPLICIT NONE 22 22 REAL(rstd) :: X=1 -
codes/icosagcm/trunk/src/etat0.f90
r170 r186 1 1 MODULE etat0_mod 2 2 CHARACTER(len=255),SAVE :: etat0_type 3 !$OMP THREADPRIVATE(etat0_type) 3 4 4 5 CONTAINS … … 26 27 TYPE(t_field),POINTER :: f_u(:) 27 28 TYPE(t_field),POINTER :: f_q(:) 29 28 30 REAL(rstd),POINTER :: ps(:), mass(:,:) 29 31 LOGICAL :: init_mass … … 75 77 76 78 IF(init_mass) THEN ! initialize mass distribution using ps 77 !$OMP BARRIER79 ! !$OMP BARRIER 78 80 DO ind=1,ndomain 81 IF (.NOT. assigned_domain(ind)) CYCLE 79 82 CALL swap_dimensions(ind) 80 83 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/etat0_academic.f90
r19 r186 10 10 USE kinetic_mod 11 11 IMPLICIT NONE 12 TYPE(t_field),POINTER :: f_ps(:)13 TYPE(t_field),POINTER :: f_phis(:)14 TYPE(t_field),POINTER :: f_theta_rhodz(:)15 TYPE(t_field),POINTER :: f_u(:)16 TYPE(t_field),POINTER :: f_q(:)17 TYPE(t_field),POINTER :: f_Ki(:)18 TYPE(t_field),POINTER :: f_temp(:)12 TYPE(t_field),POINTER,SAVE :: f_ps(:) 13 TYPE(t_field),POINTER,SAVE :: f_phis(:) 14 TYPE(t_field),POINTER,SAVE :: f_theta_rhodz(:) 15 TYPE(t_field),POINTER,SAVE :: f_u(:) 16 TYPE(t_field),POINTER,SAVE :: f_q(:) 17 TYPE(t_field),POINTER,SAVE :: f_Ki(:) 18 TYPE(t_field),POINTER,SAVE :: f_temp(:) 19 19 20 20 REAL(rstd),POINTER :: Ki(:,:) … … 58 58 59 59 DO ind=1,ndomain 60 IF (.NOT. assigned_domain(ind)) CYCLE 60 61 CALL swap_dimensions(ind) 61 62 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/etat0_dcmip1.f90
r72 r186 4 4 5 5 REAL(rstd), SAVE :: h0=1. 6 !$OMP THREADPRIVATE(h0) 6 7 REAL(rstd), SAVE :: lon0=3*pi/2 8 !$OMP THREADPRIVATE(lon0) 7 9 REAL(rstd), SAVE :: lat0=0.0 10 !$OMP THREADPRIVATE(lat0) 8 11 REAL(rstd), SAVE :: alpha=0.0 12 !$OMP THREADPRIVATE(alpha) 9 13 REAL(rstd), SAVE :: R0 14 !$OMP THREADPRIVATE(R0) 10 15 REAL(rstd), SAVE :: lat1=0. 16 !$OMP THREADPRIVATE(lat1) 11 17 REAL(rstd), SAVE :: lat2=0. 18 !$OMP THREADPRIVATE(lat2) 12 19 REAL(rstd), SAVE :: lon1=pi/6 20 !$OMP THREADPRIVATE(lon1) 13 21 REAL(rstd), SAVE :: lon2=-pi/6 22 !$OMP THREADPRIVATE(lon2) 14 23 REAL(rstd), SAVE :: latc1=0. 24 !$OMP THREADPRIVATE(latc1) 15 25 REAL(rstd), SAVE :: latc2=0. 26 !$OMP THREADPRIVATE(latc2) 16 27 REAL(rstd), SAVE :: lonc1=5*pi/6 28 !$OMP THREADPRIVATE(lonc1) 17 29 REAL(rstd), SAVE :: lonc2=7*pi/6 30 !$OMP THREADPRIVATE(lonc2) 18 31 REAL(rstd), SAVE :: zt=1000.0 32 !$OMP THREADPRIVATE(zt) 19 33 REAL(rstd), SAVE :: rt 34 !$OMP THREADPRIVATE(rt) 20 35 REAL(rstd), SAVE :: zc=5000.0 36 !$OMP THREADPRIVATE(zc) 21 37 22 38 PUBLIC etat0 … … 47 63 48 64 DO ind=1,ndomain 65 IF (.NOT. assigned_domain(ind)) CYCLE 49 66 CALL swap_dimensions(ind) 50 67 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/etat0_dcmip2.f90
r117 r186 46 46 47 47 DO ind=1,ndomain 48 IF (.NOT. assigned_domain(ind)) CYCLE 48 49 CALL swap_dimensions(ind) 49 50 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/etat0_dcmip3.f90
r55 r186 34 34 35 35 DO ind=1,ndomain 36 IF (.NOT. assigned_domain(ind)) CYCLE 36 37 CALL swap_dimensions(ind) 37 38 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/etat0_dcmip4.f90
r138 r186 12 12 REAL(rstd),PARAMETER :: up0=1 13 13 REAL(rstd) :: latw=2*Pi/9 14 !$OMP THREADPRIVATE(latw) 14 15 REAL(rstd) :: pw=34000 16 !$OMP THREADPRIVATE(pw) 15 17 REAL(rstd) :: q0=0.021 18 !$OMP THREADPRIVATE(q0) 16 19 REAL(rstd) :: lonc 20 !$OMP THREADPRIVATE(lonc) 17 21 REAL(rstd) :: latc 22 !$OMP THREADPRIVATE(latc) 23 24 INTEGER,SAVE :: testcase 25 !$OMP THREADPRIVATE(testcase) 26 18 27 PUBLIC etat0 19 28 CONTAINS … … 36 45 REAL(rstd),POINTER :: q(:,:,:) 37 46 INTEGER :: ind 47 48 testcase=1 49 CALL getin("dcmip4_testcase",testcase) 38 50 39 51 DO ind=1,ndomain 52 IF (.NOT. assigned_domain(ind)) CYCLE 40 53 CALL swap_dimensions(ind) 41 54 CALL swap_geometry(ind) … … 81 94 REAL(rstd) :: lonx,latx 82 95 REAL(rstd) :: dthetaodeta_ave, dthetaodeta, dthetaodlat, duodeta, K, r 83 INTEGER :: testcase84 96 85 97 lonc=Pi/9 86 98 latc=2*Pi/9 87 88 testcase=189 CALL getin("dcmip4_testcase",testcase)90 91 99 92 100 DO l=1,llm -
codes/icosagcm/trunk/src/etat0_dcmip5.f90
r114 r186 16 16 REAL(rstd),PARAMETER :: epsilon=1e-25 17 17 REAL(rstd),PARAMETER :: Rd=287 18 REAL(rstd) :: lonc 19 REAL(rstd) :: latc 18 REAL(rstd),SAVE :: lonc 19 !$OMP THREADPRIVATE(lonc) 20 REAL(rstd),SAVE :: latc 21 !$OMP THREADPRIVATE(latc) 20 22 TYPE(t_field),POINTER :: f_out_i(:) 21 REAL(rstd),POINTER :: out_i(:,:) 23 REAL(rstd),POINTER,SAVE :: out_i(:,:) 24 !$OMP THREADPRIVATE(out_i) 22 25 23 26 PUBLIC etat0 … … 45 48 46 49 DO ind=1,ndomain 50 IF (.NOT. assigned_domain(ind)) CYCLE 47 51 CALL swap_dimensions(ind) 48 52 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/etat0_heldsz.f90
r170 r186 8 8 TYPE(t_field),POINTER :: f_clat(:) ! FIXME, duplication 9 9 10 REAL(rstd),ALLOCATABLE :: knewt_t(:),kfrict(:)11 10 REAL(rstd),ALLOCATABLE,SAVE :: knewt_t(:),kfrict(:) 11 !$OMP THREADPRIVATE(knewt_t,kfrict) 12 12 LOGICAL, SAVE :: done=.FALSE. 13 14 REAL(rstd) :: teta0,ttp,delt_y,delt_z,eps 15 REAL(rstd) :: knewt_g, k_f,k_c_a,k_c_s 13 !$OMP THREADPRIVATE(done) 14 15 REAL(rstd),SAVE :: teta0,ttp,delt_y,delt_z,eps 16 !$OMP THREADPRIVATE(teta0,ttp,delt_y,delt_z,eps) 17 18 REAL(rstd),SAVE :: knewt_g, k_f,k_c_a,k_c_s 19 !$OMP THREADPRIVATE(knewt_g, k_f,k_c_a,k_c_s) 16 20 17 21 PUBLIC :: etat0, held_suarez … … 69 73 CALL Init_Teq 70 74 DO ind=1,ndomain 75 IF (.NOT. assigned_domain(ind)) CYCLE 71 76 CALL swap_dimensions(ind) 72 77 CALL swap_geometry(ind) … … 139 144 140 145 DO ind=1,ndomain 146 IF (.NOT. assigned_domain(ind)) CYCLE 141 147 CALL swap_dimensions(ind) 142 148 CALL swap_geometry(ind) … … 230 236 231 237 DO ind=1,ndomain 238 IF (.NOT. assigned_domain(ind)) CYCLE 232 239 CALL swap_dimensions(ind) 233 240 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/etat0_jablonowsky06.f90
r19 r186 23 23 USE vorticity_mod 24 24 IMPLICIT NONE 25 TYPE(t_field),POINTER :: f_ps(:)26 TYPE(t_field),POINTER :: f_phis(:)27 TYPE(t_field),POINTER :: f_theta_rhodz(:)28 TYPE(t_field),POINTER :: f_u(:)29 TYPE(t_field),POINTER :: f_Ki(:)30 TYPE(t_field),POINTER :: f_temp(:)31 TYPE(t_field),POINTER :: f_p(:)32 TYPE(t_field),POINTER :: f_pks(:)33 TYPE(t_field),POINTER :: f_pk(:)34 TYPE(t_field),POINTER :: f_phi(:)35 TYPE(t_field),POINTER :: f_vort(:)36 TYPE(t_field),POINTER :: f_q(:)25 TYPE(t_field),POINTER,SAVE :: f_ps(:) 26 TYPE(t_field),POINTER,SAVE :: f_phis(:) 27 TYPE(t_field),POINTER,SAVE :: f_theta_rhodz(:) 28 TYPE(t_field),POINTER,SAVE :: f_u(:) 29 TYPE(t_field),POINTER,SAVE :: f_Ki(:) 30 TYPE(t_field),POINTER,SAVE :: f_temp(:) 31 TYPE(t_field),POINTER,SAVE :: f_p(:) 32 TYPE(t_field),POINTER,SAVE :: f_pks(:) 33 TYPE(t_field),POINTER,SAVE :: f_pk(:) 34 TYPE(t_field),POINTER,SAVE :: f_phi(:) 35 TYPE(t_field),POINTER,SAVE :: f_vort(:) 36 TYPE(t_field),POINTER,SAVE :: f_q(:) 37 37 38 38 REAL(rstd),POINTER :: Ki(:,:) … … 90 90 91 91 DO ind=1,ndomain 92 IF (.NOT. assigned_domain(ind)) CYCLE 92 93 CALL swap_dimensions(ind) 93 94 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/etat0_williamson.f90
r179 r186 23 23 24 24 DO ind=1,ndomain 25 IF (.NOT. assigned_domain(ind)) CYCLE 25 26 CALL swap_dimensions(ind) 26 27 CALL swap_geometry(ind) … … 61 62 62 63 DO ind=1,ndomain 64 IF (.NOT. assigned_domain(ind)) CYCLE 63 65 CALL swap_dimensions(ind) 64 66 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/exner.f90
r151 r186 1 1 MODULE exner_mod 2 2 3 INTEGER :: caldyn_exner 3 INTEGER,SAVE :: caldyn_exner 4 !$OMP THREADPRIVATE(caldyn_exner) 5 4 6 INTEGER, PARAMETER :: lmdz=3, direct=4 5 7 … … 21 23 22 24 DO ind=1,ndomain 25 IF (.NOT. assigned_domain(ind)) CYCLE 23 26 CALL swap_dimensions(ind) 24 27 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/field.f90
r159 r186 57 57 INTEGER :: ii_size,jj_size 58 58 59 !$OMP BARRIER 60 !$OMP MASTER 59 61 ALLOCATE(field(ndomain)) 62 !$OMP END MASTER 63 !$OMP BARRIER 60 64 61 65 DO ind=1,ndomain 66 IF (.NOT. assigned_domain(ind)) CYCLE 67 62 68 IF(PRESENT(name)) THEN 63 69 field(ind)%name = name … … 106 112 107 113 ENDDO 114 !$OMP BARRIER 108 115 109 116 END SUBROUTINE allocate_field -
codes/icosagcm/trunk/src/geometry.f90
r155 r186 25 25 END TYPE t_geometry 26 26 27 TYPE(t_geometry),TARGET :: geom 27 TYPE(t_geometry),SAVE,TARGET :: geom 28 28 29 29 30 REAL(rstd),POINTER :: Ai(:) ! area of a cell 31 !$OMP THREADPRIVATE(Ai) 30 32 REAL(rstd),POINTER :: xyz_i(:,:) ! coordinate of the center of the cell (voronoi) 33 !$OMP THREADPRIVATE(xyz_i) 31 34 REAL(rstd),POINTER :: centroid(:,:) ! coordinate of the centroid of the cell 35 !$OMP THREADPRIVATE(centroid) 32 36 REAL(rstd),POINTER :: xyz_e(:,:) ! coordinate of a wind point on the cell on a edge 37 !$OMP THREADPRIVATE(xyz_e) 33 38 REAL(rstd),POINTER :: ep_e(:,:) ! perpendicular unit vector of a edge (outsider) 39 !$OMP THREADPRIVATE(ep_e) 34 40 REAL(rstd),POINTER :: et_e(:,:) ! tangeantial unit vector of a edge 41 !$OMP THREADPRIVATE(et_e) 35 42 REAL(rstd),POINTER :: elon_i(:,:) ! unit longitude vector on the center 43 !$OMP THREADPRIVATE(elon_i) 36 44 REAL(rstd),POINTER :: elat_i(:,:) ! unit latitude vector on the center 45 !$OMP THREADPRIVATE(elat_i) 37 46 REAL(rstd),POINTER :: elon_e(:,:) ! unit longitude vector on a wind point 47 !$OMP THREADPRIVATE(elon_e) 38 48 REAL(rstd),POINTER :: elat_e(:,:) ! unit latitude vector on a wind point 49 !$OMP THREADPRIVATE(elat_e) 39 50 REAL(rstd),POINTER :: xyz_v(:,:) ! coordinate of a vertex (center of the dual mesh) 51 !$OMP THREADPRIVATE(xyz_v) 40 52 REAL(rstd),POINTER :: Av(:) ! area of dual mesk cell 53 !$OMP THREADPRIVATE(Av) 41 54 REAL(rstd),POINTER :: de(:) ! distance from a neighbour == lenght of an edge of the dual mesh 55 !$OMP THREADPRIVATE(de) 42 56 REAL(rstd),POINTER :: le(:) ! lenght of a edge 57 !$OMP THREADPRIVATE(le) 43 58 REAL(rstd),POINTER :: Riv(:,:) ! weight 59 !$OMP THREADPRIVATE(Riv) 44 60 REAL(rstd),POINTER :: Riv2(:,:) ! weight 61 !$OMP THREADPRIVATE(Riv2) 45 62 INTEGER,POINTER :: ne(:,:) ! convention for the way on the normal wind on an edge 63 !$OMP THREADPRIVATE(ne) 46 64 REAL(rstd),POINTER :: Wee(:,:,:) ! weight 65 !$OMP THREADPRIVATE(Wee) 47 66 REAL(rstd),POINTER :: bi(:) ! orographie 67 !$OMP THREADPRIVATE(bi) 48 68 REAL(rstd),POINTER :: fv(:) ! coriolis (evaluted on a vertex) 69 !$OMP THREADPRIVATE(fv) 49 70 50 71 … … 91 112 IMPLICIT NONE 92 113 INTEGER,INTENT(IN) :: ind 93 ! $OMP MASTER114 !!$OMP MASTER 94 115 Ai=geom%Ai(ind) 95 116 xyz_i=geom%xyz_i(ind) … … 112 133 bi=geom%bi(ind) 113 134 fv=geom%fv(ind) 114 ! $OMP END MASTER115 ! $OMP BARRIER135 !!$OMP END MASTER 136 !!$OMP BARRIER 116 137 END SUBROUTINE swap_geometry 117 138 … … 127 148 REAL(rstd) :: vect(3,6) 128 149 REAL(rstd) :: centr(3) 129 130 150 INTEGER :: ind,i,j,n,k 131 132 CALL transfert_request(geom%xyz_i,req_i1) 151 TYPE(t_message),SAVE :: message 152 LOGICAL, SAVE :: first=.TRUE. 153 !$OMP THREADPRIVATE(first) 154 155 IF (first) THEN 156 CALL init_message(geom%xyz_i, req_i1 ,message) 157 first=.FALSE. 158 ENDIF 159 160 CALL transfert_message(geom%xyz_i,message) 161 162 163 ! CALL transfert_request(geom%xyz_i,req_i1) 133 164 134 165 DO ind=1,ndomain 166 IF (.NOT. assigned_domain(ind)) CYCLE 135 167 CALL swap_dimensions(ind) 136 168 CALL swap_geometry(ind) … … 157 189 USE transfert_mod 158 190 USE vector 159 USE ioipsl191 USE getin_mod 160 192 IMPLICIT NONE 161 193 INTEGER :: nb_it=0 … … 172 204 173 205 DO ind=1,ndomain 206 IF (.NOT. assigned_domain(ind)) CYCLE 174 207 d=>domain(ind) 175 208 CALL swap_dimensions(ind) … … 186 219 187 220 DO ind=1,ndomain 221 IF (.NOT. assigned_domain(ind)) CYCLE 188 222 d=>domain(ind) 189 223 CALL swap_dimensions(ind) … … 215 249 sum=0 216 250 DO ind=1,ndomain 251 IF (.NOT. assigned_domain(ind)) CYCLE 217 252 CALL swap_dimensions(ind) 218 253 CALL swap_geometry(ind) … … 253 288 USE dimensions 254 289 USE transfert_mod 255 USE ioipsl290 USE getin_mod 256 291 IMPLICIT NONE 257 292 … … 273 308 274 309 DO ind=1,ndomain 310 IF (.NOT. assigned_domain(ind)) CYCLE 275 311 d=>domain(ind) 276 312 CALL swap_dimensions(ind) … … 489 525 490 526 ENDDO 527 491 528 CALL transfert_request(geom%Ai,req_i1) 492 529 CALL transfert_request(geom%centroid,req_i1) 530 493 531 CALL surf_triangle(d%xyz(:,ii_begin,jj_begin),d%xyz(:,ii_begin,jj_end),d%xyz(:,ii_end,jj_begin),S) 494 532 -
codes/icosagcm/trunk/src/geopotential_mod.f90
r121 r186 20 20 21 21 DO ind=1,ndomain 22 IF (.NOT. assigned_domain(ind)) CYCLE 22 23 CALL swap_dimensions(ind) 23 24 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/guided_mod.f90
r98 r186 2 2 3 3 CHARACTER(LEN=255),SAVE :: guided_type 4 !$OMP THREADPRIVATE(guided_type) 4 5 5 6 CONTAINS -
codes/icosagcm/trunk/src/guided_ncar_mod.f90
r111 r186 4 4 5 5 INTEGER,SAVE :: case_wind 6 !$OMP THREADPRIVATE(case_wind) 6 7 7 8 REAL(rstd), PARAMETER :: alpha=0.0 ! tilt of solid-body rotation … … 46 47 47 48 DO ind = 1 , ndomain 49 IF (.NOT. assigned_domain(ind)) CYCLE 48 50 CALL swap_dimensions(ind) 49 51 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/icosa_gcm.f90
r171 r186 12 12 USE xios_mod 13 13 USE write_field 14 USE getin_mod 14 15 IMPLICIT NONE 15 16 16 TYPE(t_field),POINTER :: sum_ne(:)17 TYPE(t_field),POINTER :: sum_ne_glo(:)17 TYPE(t_field),POINTER,SAVE :: sum_ne(:) 18 TYPE(t_field),POINTER,SAVE :: sum_ne_glo(:) 18 19 REAL(rstd),POINTER :: pt_sum_ne(:) 19 20 … … 24 25 25 26 CALL init_mpipara 27 CALL trace_off 26 28 CALL xios_init 27 29 CALL init_earth_const … … 34 36 CALL init_trace 35 37 38 !$OMP PARALLEL 39 CALL compute_geometry 36 40 37 CALL compute_geometry38 41 CALL init_disvert 39 42 CALL init_vertical_interp … … 41 44 CALL allocate_field(sum_ne,field_T,type_real) 42 45 43 46 !$OMP BARRIER 47 !$OMP MASTER 44 48 DO ind=1,ndomain 45 49 … … 60 64 ENDDO 61 65 ENDDO 62 63 64 66 65 67 IF (is_mpi_root) PRINT *," Diff surf",1-tot_sum/(4*Pi*radius*radius) 68 !$OMP END MASTER 66 69 70 71 CALL WriteField("Ai",geom%Ai) 67 72 68 CALL WriteField("Ai",geom%Ai)69 70 73 IF (is_mpi_root) CALL write_apbp 71 74 CALL init_time … … 74 77 CALL init_timeloop 75 78 76 !$OMP PARALLEL 79 77 80 CALL timeloop 81 78 82 !$OMP END PARALLEL 79 83 -
codes/icosagcm/trunk/src/icosa_mod.f90
r82 r186 2 2 3 3 USE genmod 4 USE ioipsl, ONLY : getin 4 ! USE ioipsl, ONLY : getin 5 USE getin_mod, ONLY : getin 5 6 USE grid_param 6 7 USE metric -
codes/icosagcm/trunk/src/kinetic.f90
r146 r186 18 18 19 19 DO ind=1,ndomain 20 IF (.NOT. assigned_domain(ind)) CYCLE 20 21 CALL swap_dimensions(ind) 21 22 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/mpi_mod.F90
r172 r186 12 12 INTEGER :: MPI_STATUS_SIZE 13 13 INTEGER :: MPI_SUM 14 INTEGER :: MPI_THREAD_SINGLE, MPI_THREAD_FUNNELED, 15 INTEGER :: MPI_THREAD_SERIALIZED, MPI_THREAD_MULTIPLE 16 14 17 INTEGER,PARAMETER :: MPI_ADDRESS_KIND=KIND(INTEGER) 15 18 #endif … … 21 24 22 25 SUBROUTINE MPI_INIT 26 PRINT *, 'Compiled without MPI' 27 END 28 29 SUBROUTINE MPI_INIT_THREAD 23 30 PRINT *, 'Compiled without MPI' 24 31 END … … 69 76 END 70 77 78 SUBROUTINE MPI_FREE_MEM 79 END 80 71 81 #endif -
codes/icosagcm/trunk/src/mpipara.F90
r171 r186 3 3 INTEGER,SAVE :: mpi_rank 4 4 INTEGER,SAVE :: mpi_size 5 INTEGER,SAVE :: mpi_threading_mode 5 6 6 7 INTEGER,SAVE :: comm_icosa … … 13 14 END INTERFACE allocate_mpi_buffer 14 15 16 INTERFACE free_mpi_buffer 17 MODULE PROCEDURE free_mpi_buffer_r2, free_mpi_buffer_r3, free_mpi_buffer_r4 18 END INTERFACE free_mpi_buffer 19 15 20 CONTAINS 16 21 17 22 SUBROUTINE init_mpipara 18 23 USE mpi_mod 24 USE getin_mod 19 25 #ifdef CPP_USING_XIOS 20 26 USE xios 21 27 #endif 22 28 IMPLICIT NONE 29 CHARACTER(LEN=256) :: required_mode_str 30 INTEGER :: required_mode 23 31 24 32 using_mpi=.FALSE. … … 28 36 29 37 IF (using_mpi) THEN 30 CALL MPI_INIT(ierr) 31 38 39 required_mode_str='multiple' 40 CALL getin('mpi_threading_mode',required_mode_str) 41 42 SELECT CASE(TRIM(required_mode_str)) 43 CASE ('single') 44 required_mode=MPI_THREAD_SINGLE 45 CASE ('funneled') 46 required_mode=MPI_THREAD_FUNNELED 47 CASE ('serialized') 48 required_mode=MPI_THREAD_SERIALIZED 49 CASE ('multiple') 50 required_mode=MPI_THREAD_MULTIPLE 51 CASE DEFAULT 52 PRINT*,'Bad selector for variable mpi_threading_mode : <', TRIM(required_mode_str), & 53 '> => options are <single>, <funneled>, <serialized>, <multiple>' 54 STOP 55 END SELECT 56 57 58 IF (required_mode==MPI_THREAD_SINGLE) PRINT*,'MPI_INIT_THREAD : MPI_SINGLE_THREAD required' 59 IF (required_mode==MPI_THREAD_FUNNELED) PRINT*,'MPI_INIT_THREAD : MPI_THREAD_FUNNELED required' 60 IF (required_mode==MPI_THREAD_SERIALIZED) PRINT*,'MPI_INIT_THREAD : MPI_THREAD_SERIALIZED required' 61 IF (required_mode==MPI_THREAD_MULTIPLE) PRINT*,'MPI_INIT_THREAD : MPI_THREAD_MULTIPLE required' 62 63 64 CALL MPI_INIT_THREAD(MPI_THREAD_MULTIPLE,mpi_threading_mode,ierr) 65 66 IF (mpi_threading_mode==MPI_THREAD_SINGLE) PRINT*,'MPI_INIT_THREAD : MPI_SINGLE_THREAD provided' 67 IF (mpi_threading_mode==MPI_THREAD_FUNNELED) PRINT*,'MPI_INIT_THREAD : MPI_THREAD_FUNNELED provided' 68 IF (mpi_threading_mode==MPI_THREAD_SERIALIZED) PRINT*,'MPI_INIT_THREAD : MPI_THREAD_SERIALIZED provided' 69 IF (mpi_threading_mode==MPI_THREAD_MULTIPLE) PRINT*,'MPI_INIT_THREAD : MPI_THREAD_MULTIPLE provided' 70 71 IF (mpi_threading_mode > required_mode) mpi_threading_mode=required_mode 72 73 IF (mpi_threading_mode==MPI_THREAD_SINGLE) PRINT*,'MPI_INIT_THREAD : MPI_SINGLE_THREAD used' 74 IF (mpi_threading_mode==MPI_THREAD_FUNNELED) PRINT*,'MPI_INIT_THREAD : MPI_THREAD_FUNNELED used' 75 IF (mpi_threading_mode==MPI_THREAD_SERIALIZED) PRINT*,'MPI_INIT_THREAD : MPI_THREAD_SERIALIZED used' 76 IF (mpi_threading_mode==MPI_THREAD_MULTIPLE) PRINT*,'MPI_INIT_THREAD : MPI_THREAD_MULTIPLE used' 77 32 78 #ifdef CPP_USING_XIOS 33 79 CALL xios_initialize("icosagcm",return_comm=comm_icosa) … … 79 125 CALL C_F_POINTER(base_ptr, buffer, (/ length /)) 80 126 81 END SUBROUTINE allocate_mpi_buffer_r2 127 END SUBROUTINE allocate_mpi_buffer_r2 128 129 SUBROUTINE free_mpi_buffer_r2(buffer) 130 USE ISO_C_BINDING 131 USE mpi_mod 132 USE prec 133 IMPLICIT NONE 134 REAL(rstd), POINTER :: buffer(:) 135 136 CALL MPI_FREE_MEM(buffer,ierr) 137 138 END SUBROUTINE free_mpi_buffer_r2 82 139 83 140 SUBROUTINE allocate_mpi_buffer_r3(buffer,length,dim3) … … 100 157 CALL C_F_POINTER(base_ptr, buffer, (/ length,dim3 /)) 101 158 102 END SUBROUTINE allocate_mpi_buffer_r3 159 END SUBROUTINE allocate_mpi_buffer_r3 160 161 SUBROUTINE free_mpi_buffer_r3(buffer) 162 USE ISO_C_BINDING 163 USE mpi_mod 164 USE prec 165 IMPLICIT NONE 166 REAL(rstd), POINTER :: buffer(:,:) 167 168 CALL MPI_FREE_MEM(buffer,ierr) 169 170 END SUBROUTINE free_mpi_buffer_r3 103 171 104 172 SUBROUTINE allocate_mpi_buffer_r4(buffer,length,dim3,dim4) … … 123 191 124 192 END SUBROUTINE allocate_mpi_buffer_r4 193 194 SUBROUTINE free_mpi_buffer_r4(buffer) 195 USE ISO_C_BINDING 196 USE mpi_mod 197 USE prec 198 IMPLICIT NONE 199 REAL(rstd), POINTER :: buffer(:,:,:) 200 201 CALL MPI_FREE_MEM(buffer,ierr) 202 203 END SUBROUTINE free_mpi_buffer_r4 125 204 126 205 END MODULE mpipara -
codes/icosagcm/trunk/src/omega.f90
r96 r186 13 13 REAL(rstd),POINTER :: ps(:), u(:,:), om(:,:) 14 14 DO ind=1,ndomain 15 IF (.NOT. assigned_domain(ind)) CYCLE 15 16 CALL swap_dimensions(ind) 16 17 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/omp_para.F90
r151 r186 18 18 LOGICAL,SAVE :: using_openmp 19 19 20 LOGICAL,PARAMETER :: omp_by_domain=.TRUE. 21 20 22 CONTAINS 21 23 … … 23 25 SUBROUTINE init_omp_para 24 26 USE grid_param 25 #ifdef CPP_USING_O PENMP27 #ifdef CPP_USING_OMP 26 28 USE omp_lib 27 29 #endif … … 29 31 INTEGER :: ll_nb,i 30 32 31 #ifdef CPP_USING_O PENMP33 #ifdef CPP_USING_OMP 32 34 using_openmp=.TRUE. 33 35 #else … … 39 41 40 42 !$OMP MASTER 41 #ifdef CPP_USING_O PENMP43 #ifdef CPP_USING_OMP 42 44 omp_size=OMP_GET_NUM_THREADS() 43 45 #endif 44 46 !$OMP END MASTER 45 47 !$OMP BARRIER 46 #ifdef CPP_USING_O PENMP48 #ifdef CPP_USING_OMP 47 49 omp_rank=OMP_GET_THREAD_NUM() 48 50 #endif 49 omp_first=.FALSE. 50 omp_last=.FALSE. 51 omp_master=.FALSE. 51 52 IF (omp_by_domain) THEN 53 omp_first=.TRUE. 54 omp_last=.TRUE. 55 IF (omp_rank==0) THEN 56 omp_master=.TRUE. 57 ELSE 58 omp_master=.FALSE. 59 ENDIF 60 61 ll_begin=1 62 ll_beginp1=2 63 ll_end=llm 64 ll_endm1=llm-1 65 ll_endp1=llm+1 66 67 ELSE 52 68 53 IF (omp_rank==0) THEN 54 omp_first=.TRUE. 55 omp_master=.TRUE. 69 omp_first=.FALSE. 70 omp_last=.FALSE. 71 omp_master=.FALSE. 72 73 IF (omp_rank==0) THEN 74 omp_first=.TRUE. 75 omp_master=.TRUE. 76 ENDIF 77 78 IF (omp_rank==omp_size-1) omp_last=.TRUE. 79 80 ll_end=0 81 DO i=0,omp_rank 82 ll_begin=ll_end+1 83 ll_nb=llm/omp_size 84 IF (MOD(llm,omp_size)>i) ll_nb=ll_nb+1 85 ll_end=ll_begin+ll_nb-1 86 ENDDO 87 88 ll_beginp1=ll_begin 89 ll_endp1=ll_end 90 ll_endm1=ll_end 91 92 IF (omp_first) ll_beginp1=ll_begin+1 93 IF (omp_last) ll_endp1=ll_endp1+1 94 IF (omp_last) ll_endm1=ll_endm1-1 95 56 96 ENDIF 57 58 IF (omp_rank==omp_size-1) omp_last=.TRUE.59 60 ll_end=061 DO i=0,omp_rank62 ll_begin=ll_end+163 ll_nb=llm/omp_size64 IF (MOD(llm,omp_size)>i) ll_nb=ll_nb+165 ll_end=ll_begin+ll_nb-166 ENDDO67 68 ll_beginp1=ll_begin69 ll_endp1=ll_end70 ll_endm1=ll_end71 72 IF (omp_first) ll_beginp1=ll_begin+173 IF (omp_last) ll_endp1=ll_endp1+174 IF (omp_last) ll_endm1=ll_endm1-175 76 97 !$OMP END PARALLEL 77 98 … … 91 112 END SUBROUTINE init_omp_para 92 113 93 114 FUNCTION omp_in_parallel() 115 #ifdef CPP_USING_OMP 116 USE omp_lib, ONLY : omp_in_parallel_=>omp_in_parallel 117 #endif 118 IMPLICIT NONE 119 LOGICAL :: omp_in_parallel 120 121 #ifdef CPP_USING_OMP 122 omp_in_parallel=omp_in_parallel_() 123 #else 124 omp_in_parallel=.FALSE. 125 #endif 126 127 END FUNCTION omp_in_parallel 128 94 129 END MODULE omp_para 95 130 -
codes/icosagcm/trunk/src/output_field.f90
r171 r186 2 2 USE genmod 3 3 PRIVATE 4 LOGICAL,SAVE :: xios_output 4 LOGICAL,SAVE :: xios_output 5 !$OMP THREADPRIVATE(xios_output) 6 LOGICAL,SAVE :: enable_io 7 !$OMP THREADPRIVATE(enable_io) 5 8 6 PUBLIC xios_output,output_field_init,output_field,output_field_finalize9 PUBLIC enable_io,xios_output,output_field_init,output_field,output_field_finalize 7 10 8 11 CONTAINS 9 12 10 13 SUBROUTINE output_field_init 11 USE IOIPSL14 USE getin_mod 12 15 USE xios_mod 13 16 USE write_field 14 17 IMPLICIT NONE 18 19 enable_io=.TRUE. 20 CALL getin('enable_io',enable_io) 15 21 16 22 IF (using_xios) THEN -
codes/icosagcm/trunk/src/phyparam.F
r163 r186 1 MODULE PHY2 USE dimphys_mod3 LOGICAL:: firstcall,lastcall4 contains5 6 SUBROUTINE phyparam_lmd(it,ngrid,nlayer,nq,7 s ptimestep,lati,long,rjourvrai,gmtime,8 s pplev,pplay,pphi,pphis,9 s pu,pv,pt,pq,10 s pdu,pdv,pdt,pdq,pdpsrf)11 12 USE ICOSA13 USE dimphys_mod14 USE RADIATION15 USE SURFACE_PROCESS16 c17 IMPLICIT NONE18 c=======================================================================19 c20 c subject:21 c --------22 c23 c Organisation of the physical parametrisations of the LMD24 c 20 parameters GCM for planetary atmospheres.25 c It includes:26 c raditive transfer (long and shortwave) for CO2 and dust.27 c vertical turbulent mixing28 c convective adjsutment29 c30 c author: Frederic Hourdin 15 / 10 /9331 c -------32 c33 c arguments:34 c ----------35 c36 c input:37 c ------38 c39 c ngrid Size of the horizontal grid.40 c All internal loops are performed on that grid.41 c nlayer Number of vertical layers.42 c nq Number of advected fields43 c firstcall True at the first call44 c lastcall True at the last call45 c rjourvrai Number of days counted from the North. Spring46 c equinoxe.47 c gmtime hour (s)48 c ptimestep timestep (s)49 c pplay(ngrid,nlayer+1) Pressure at the middle of the layers (Pa)50 c pplev(ngrid,nlayer+1) intermediate pressure levels (pa)51 c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)52 c pu(ngrid,nlayer) u component of the wind (ms-1)53 c pv(ngrid,nlayer) v component of the wind (ms-1)54 c pt(ngrid,nlayer) Temperature (K)55 c pq(ngrid,nlayer,nq) Advected fields56 c pudyn(ngrid,nlayer) \57 c pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the58 c ptdyn(ngrid,nlayer) / corresponding variables59 c pqdyn(ngrid,nlayer,nq) /60 c pw(ngrid,?) vertical velocity61 c62 c output:63 c -------64 c65 c pdu(ngrid,nlayer) \66 c pdv(ngrid,nlayer) \ Temporal derivative of the corresponding67 c pdt(ngrid,nlayer) / variables due to physical processes.68 c pdq(ngrid,nlayer) /69 c pdpsrf(ngrid) /70 c71 c=======================================================================72 c73 !c-----------------------------------------------------------------------74 !c75 !c 0. Declarations :76 !c ------------------77 !c Arguments :78 !c -----------79 80 !c inputs:81 !c -------82 INTEGER ngrid,nlayer,nq,it,ij,i83 REAL ptimestep84 real zdtime85 REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)86 REAL pphi(ngrid,nlayer)87 REAL pphis(ngrid)88 REAL pu(ngrid,nlayer),pv(ngrid,nlayer)89 REAL pt(ngrid,nlayer),pq(ngrid,nlayer,nq)90 REAL pdu(ngrid,nlayer),pdv(ngrid,nlayer)91 92 !c dynamial tendencies93 REAL pdtdyn(ngrid,nlayer),pdqdyn(ngrid,nlayer,nq)94 REAL pdudyn(ngrid,nlayer),pdvdyn(ngrid,nlayer)95 REAL pw(ngrid,nlayer)96 97 !c Time98 real rjourvrai99 REAL gmtime100 101 !c outputs:102 !c --------103 104 !c physical tendencies105 REAL pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq)106 REAL pdpsrf(ngrid)107 108 109 !c Local variables :110 !c -----------------111 112 INTEGER j,l,ig,ierr,aslun,nlevel,igout,it1,it2,unit,isoil113 REAL zps_av114 REAL zday115 REAL zh(ngrid,nlayer),z1,z2116 REAL zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer)117 REAL zdvfr(ngrid,nlayer),zdufr(ngrid,nlayer)118 REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid)119 REAL zflubid(ngrid),zpmer(ngrid)120 REAL zplanck(ngrid),zpopsk(ngrid,nlayer)121 REAL zdum1(ngrid,nlayer)122 REAL zdum2(ngrid,nlayer)123 REAL zdum3(ngrid,nlayer)124 REAL ztim1,ztim2,ztim3125 REAL zls,zinsol126 REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer)127 REAL zfluxsw(ngrid),zfluxlw(ngrid)128 character*2 str2129 130 !c Local saved variables:131 !c ----------------------132 REAL(rstd)::long(ngrid),lati(ngrid)133 REAL(rstd)::mu0(ngrid),fract(ngrid),coslat(ngrid)134 REAL(rstd)::sinlon(ngrid),coslon(ngrid),sinlat(ngrid)135 REAL(rstd)::dist_sol,declin136 REAL::totarea !sarvesh137 !!!!!!!!sarvesh !!!!!!! CHECK SAVE ATTRIBUTE138 INTEGER:: icount139 real:: zday_last140 REAL:: solarcst141 142 SAVE icount,zday_last143 SAVE solarcst144 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!145 REAL stephan146 SAVE stephan147 DATA stephan/5.67e-08/148 DATA solarcst/1370./149 REAL presnivs(nlayer)150 INTEGER:: nn1,nn2151 c-----------------------------------------------------------------------152 c 1. Initialisations :153 c --------------------154 c call initial0(ngrid*nlayer*nqmx,pdq)155 156 ! nn1=(jj_begin -1)*iim+ii_begin157 ! nn2=(jj_end -1)*iim+ii_end158 159 nlevel=nlayer+1160 igout=ngrid/2+1161 162 DO j=jj_begin-offset,jj_end+offset163 DO i=ii_begin-offset,ii_end+offset164 ig=(j-1)*iim+i165 sinlat(ig) = sin(lati(ig))166 coslat(ig) = cos(lati(ig))167 sinlon(ig) = sin(long(ig))168 coslon(ig) = cos(long(ig))169 END DO170 ENDDO171 zday=rjourvrai+gmtime172 IF ( it == 0 ) then173 firstcall=.TRUE.174 ELSE175 firstcall=.FALSE.176 ENDIF177 IF ( it == ndays*day_step ) Then178 lastcall = .True.179 END IF180 181 IF(firstcall) THEN182 PRINT*,'FIRSTCALL ',ngridmx,nlayermx,nsoilmx183 zday_last=rjourvrai184 inertie=2000185 albedo=0.2186 emissiv=1.187 z0=0.1188 rnatur=1.189 q2=1.e-10190 q2l=1.e-10191 tsurf(:)=300.192 tsoil(:,:)=300.193 ! print*,tsoil(ngrid/2+1,nsoilmx/2+2)194 ! print*,'TS ',tsurf(igout),tsoil(igout,5)195 196 IF (.not.callrad) then197 DO j=jj_begin-offset,jj_end+offset198 DO i=ii_begin-offset,ii_end+offset199 ig=(j-1)*iim+i200 fluxrad(ig)=0.201 enddo202 enddo203 ENDIF204 205 IF(callsoil) THEN206 CALL soil(ngrid,nsoilmx,firstcall,inertie,207 s ptimestep,tsurf,tsoil,capcal,fluxgrd)208 ELSE209 PRINT*,'WARNING!!! Thermal conduction in the soil210 s turned off'211 DO j=jj_begin-offset,jj_end+offset212 DO i=ii_begin-offset,ii_end+offset213 ig=(j-1)*iim+i214 capcal(ig)=1.e5215 fluxgrd(ig)=0.216 ENDDO217 ENDDO218 ENDIF219 icount=0220 ENDIF221 222 IF (ngrid.NE.ngrid) THEN223 PRINT*,'STOP in inifis'224 PRINT*,'Probleme de dimenesions :'225 PRINT*,'ngrid = ',ngrid226 PRINT*,'ngrid = ',ngrid227 STOP228 ENDIF229 230 DO l=1,nlayer231 DO j=jj_begin-offset,jj_end+offset232 DO i=ii_begin-offset,ii_end+offset233 ig=(j-1)*iim+i234 pdv(ig,l)=0.235 pdu(ig,l)=0.236 pdt(ig,l)=0.237 ENDDO238 ENDDO239 ENDDO240 241 DO j=jj_begin-offset,jj_end+offset242 DO i=ii_begin-offset,ii_end+offset243 ig=(j-1)*iim+i244 pdpsrf(ig)=0.245 zflubid(ig)=0.246 zdtsrf(ig)=0.247 ENDDO248 ENDDO249 250 zps_av=0.0251 DO j=jj_begin-offset,jj_end+offset252 DO i=ii_begin-offset,ii_end+offset253 ig=(j-1)*iim+i254 zps_av=zps_av+pplev(ig,1)*Ai(ig)255 totarea=totarea+Ai(ig)256 END DO257 END DO258 zps_av=zps_av/totarea259 260 261 !print*,"maxplev",maxval(pplev(:,1)),minval(pplev(:,1))262 c-----------------------------------------------------------------------263 c calcul du geopotentiel aux niveaux intercouches264 c ponderation des altitudes au niveau des couches en dp/p265 266 DO l=1,nlayer267 DO j=jj_begin-offset,jj_end+offset268 DO i=ii_begin-offset,ii_end+offset269 ig=(j-1)*iim+i270 zzlay(ig,l)=pphi(ig,l)/g271 ENDDO272 ENDDO273 ENDDO274 275 !print*,"zzlay",maxval(zzlay(:,1)),minval(zzlay(:,1))276 277 278 DO j=jj_begin-offset,jj_end+offset279 DO i=ii_begin-offset,ii_end+offset280 ig=(j-1)*iim+i281 zzlev(ig,1)=0.282 ENDDO283 ENDDO284 285 DO l=2,nlayer286 DO j=jj_begin-offset,jj_end+offset287 DO i=ii_begin-offset,ii_end+offset288 ig=(j-1)*iim+i289 z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))290 z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))291 zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)292 ENDDO293 ENDDO294 ENDDO295 !print*,"zzlev",maxval(zzlev(:,1)),minval(zzlev(:,1))296 c-----------------------------------------------------------------------297 c Transformation de la temperature en temperature potentielle298 DO l=1,nlayer299 DO j=jj_begin-offset,jj_end+offset300 DO i=ii_begin-offset,ii_end+offset301 ig=(j-1)*iim +i302 zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**kappa303 ENDDO304 ENDDO305 ENDDO306 307 DO l=1,nlayer308 DO j=jj_begin-offset,jj_end+offset309 DO i=ii_begin-offset,ii_end+offset310 ig=(j-1)*iim+i311 zh(ig,l)=pt(ig,l)/zpopsk(ig,l)312 ENDDO313 ENDDO314 ENDDO315 !print*,"ph pot",maxval(zh(:,1)),minval(zh(:,1))316 ! go to 101317 c-----------------------------------------------------------------------318 c 2. Calcul of the radiative tendencies :319 c ---------------------------------------320 ! print*,'callrad0'321 IF(callrad) THEN322 ! print*,'callrad'323 ! WARNING !!! on calcule le ray a chaque appel324 ! IF( MOD(icount,iradia).EQ.0) THEN325 326 CALL solarlong(zday,zls)327 CALL orbite(zls,dist_sol,declin)328 !329 ! print*,'ATTENTIOn : pdeclin = 0',' L_s=',zls330 ! print*,'diurnal=',diurnal331 IF(diurnal) THEN332 if ( 1.eq.1 ) then333 ztim1=SIN(declin)334 ztim2=COS(declin)*COS(2.*pi*(zday-.5))335 ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))336 337 CALL solang(ngrid,sinlon,coslon,sinlat,coslat,338 s ztim1,ztim2,ztim3,339 s mu0,fract)340 else341 zdtime=ptimestep*float(iradia)342 ! CALL zenang(zls,gmtime,zdtime,lati,long,mu0,fract) ! FIXME343 !print*,'ZENANG '344 endif345 346 IF(lverbose) THEN347 PRINT*,'day, declin, sinlon,coslon,sinlat,coslat'348 PRINT*,zday, declin,349 s sinlon(igout),coslon(igout),350 s sinlat(igout),coslat(igout)351 ENDIF352 ELSE353 !print*,'declin,ngrid,radius',declin,ngrid,radius354 CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,radius)355 ! open(100,file="mu0.txt")356 ! write(100,*)(mu0(ij),ij=1,ngrid)357 ENDIF358 ! print*,"iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii"359 !c 2.2 Calcul of the radiative tendencies and fluxes:360 !c --------------------------------------------------361 !c 2.1.2 levels362 363 zinsol=solarcst/(dist_sol*dist_sol)364 ! print*,'iim,jjm,llm,ngrid,nlayer,ngrid,nlayer'365 ! print*,iim,jjm,llm,ngrid,nlayer,ngrid,nlayer366 ! print*,"zinsol sol_dist",zinsol,dist_sol367 ! STOP368 369 CALL sw(ngrid,nlayer,diurnal,coefvis,albedo,370 $ pplev,zps_av,371 $ mu0,fract,zinsol,372 $ zfluxsw,zdtsw,373 $ lverbose)374 375 !print*,"sw",maxval(zfluxsw),minval(zfluxsw),376 ! $ maxval(zdtsw),minval(zdtsw), " it",it377 378 ! print*,"lllllllllllllllllllllllllllllllllllllllll"379 ! print*,"pplev",maxval(pplev),minval(pplev)380 ! print*,"zps,tsurf",zps_av,maxval(tsurf),minval(tsurf)381 ! print*,"pt",maxval(pt),minval(pt)382 383 384 CALL lw(ngrid,nlayer,coefir,emissiv,385 $ pplev,zps_av,tsurf,pt,386 $ zfluxlw,zdtlw,387 $ lverbose)388 389 !print*,"lw",maxval(zfluxlw),minval(zfluxlw),390 ! $ maxval(zdtlw),minval(zdtlw)391 392 ! print*,"lw",maxval(393 394 ! print*,"after lw"395 c 2.4 total flux and tendencies:396 c ------------------------------397 398 c 2.4.1 fluxes399 400 DO j=jj_begin-offset,jj_end+offset401 DO i=ii_begin-offset,ii_end+offset402 ig=(j-1)*iim+i403 fluxrad(ig)=emissiv(ig)*zfluxlw(ig)404 $ +zfluxsw(ig)*(1.-albedo(ig))405 406 zplanck(ig)=tsurf(ig)*tsurf(ig)407 408 zplanck(ig)=emissiv(ig)*409 $ stephan*zplanck(ig)*zplanck(ig)410 411 fluxrad(ig)=fluxrad(ig)-zplanck(ig)412 ENDDO413 ENDDO414 c 2.4.2 temperature tendencies415 416 DO l=1,nlayer417 DO j=jj_begin-offset,jj_end+offset418 DO i=ii_begin-offset,ii_end+offset419 ig=(j-1)*iim+i420 dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)421 ENDDO422 ENDDO423 ENDDO424 425 426 c 2.5 Transformation of the radiative tendencies:427 c -----------------------------------------------428 429 DO l=1,nlayer430 DO j=jj_begin-offset,jj_end+offset431 DO i=ii_begin-offset,ii_end+offset432 ig=(j-1)*iim+i433 pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)434 ENDDO435 ENDDO436 ENDDO437 438 IF(lverbose) THEN439 PRINT*,'Diagnotique for the radaition'440 PRINT*,'albedo, emissiv, mu0,fract,Frad,Planck'441 PRINT*,albedo(igout),emissiv(igout),mu0(igout),442 s fract(igout),443 s fluxrad(igout),zplanck(igout)444 PRINT*,'Tlay Play Plev dT/dt SW dT/dt LW (K/day)'445 ! PRINT*,'unjours',unjours446 DO l=1,nlayer447 WRITE(*,'(3f15.5,2e15.2)') pt(igout,l),448 s pplay(igout,l),pplev(igout,l),449 s zdtsw(igout,l),zdtlw(igout,l)450 ENDDO451 ENDIF452 ENDIF !( CALL RADIATION )453 ! print*,"eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee"454 455 c-----------------------------------------------------------------------456 c 3. Vertical diffusion (turbulent mixing):457 c -----------------------------------------458 c459 IF(calldifv) THEN460 461 DO ig=1,ngrid462 zflubid(ig)=fluxrad(ig)+fluxgrd(ig)463 ENDDO464 465 CALL zerophys(ngrid*nlayer,zdum1)466 CALL zerophys(ngrid*nlayer,zdum2)467 c CALL zerophys(ngrid*nlayer,zdum3)468 do l=1,nlayer469 do ig=1,ngrid470 zdum3(ig,l)=pdt(ig,l)/zpopsk(ig,l)471 enddo472 enddo473 474 CALL vdif(ngrid,nlayer,zday,475 $ ptimestep,capcal,z0,476 $ pplay,pplev,zzlay,zzlev,477 $ pu,pv,zh,tsurf,emissiv,478 $ zdum1,zdum2,zdum3,zflubid,479 $ zdufr,zdvfr,zdhfr,zdtsrfr,q2,q2l,480 $ lverbose)481 c igout=ngrid/2+1482 c PRINT*,'zdufr zdvfr zdhfr'483 c DO l=1,nlayer484 c PRINT*,zdufr(igout,1),zdvfr(igout,l),zdhfr(igout,l)485 c ENDDO486 c CALL difv (ngrid,nlayer,487 c $ area,lati,capcal,488 c $ pplev,pphi,489 c $ pu,pv,zh,tsurf,emissiv,490 c $ zdum1,zdum2,zdum3,zflubid,491 c $ zdufr,zdvfr,zdhfr,zdtsrf,492 c $ .true.)493 c PRINT*,'zdufr zdvfr zdhfr'494 c DO l=1,nlayer495 c PRINT*,zdufr(igout,1),zdvfr(igout,l),zdhfr(igout,l)496 c ENDDO497 c STOP498 499 DO l=1,nlayer500 DO ig=1,ngrid501 pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)502 pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)503 pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)504 ENDDO505 ENDDO506 507 DO ig=1,ngrid508 zdtsrf(ig)=zdtsrf(ig)+zdtsrfr(ig)509 ENDDO510 511 ELSE512 513 DO j=jj_begin-offset,jj_end+offset514 DO i=ii_begin-offset,ii_end+offset515 ig=(j-1)*iim+i516 zdtsrf(ig)=zdtsrf(ig)+517 s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)518 ENDDO519 ENDDO520 521 ! write(66,*)"tsrf",maxval(zdtsrf(:)),minval(zdtsrf(:))522 ! write(66,*)"frd",maxval(fluxrad(:)),minval(fluxrad(:))523 ! write(66,*)"fgd",maxval(fluxgrd(:)),minval(fluxgrd(:))524 ENDIF525 c-----------------------------------------------------------------------526 c 4. Dry convective adjustment:527 c -----------------------------528 529 IF(calladj) THEN530 531 DO l=1,nlayer532 DO ig=1,ngrid533 zdum1(ig,l)=pdt(ig,l)/zpopsk(ig,l)534 ENDDO535 ENDDO536 CALL zerophys(ngrid*nlayer,zdufr)537 CALL zerophys(ngrid*nlayer,zdvfr)538 CALL zerophys(ngrid*nlayer,zdhfr)539 CALL convadj(ngrid,nlayer,ptimestep,540 $ pplay,pplev,zpopsk,541 $ pu,pv,zh,542 $ pdu,pdv,zdum1,543 $ zdufr,zdvfr,zdhfr)544 545 DO l=1,nlayer546 DO ig=1,ngrid547 pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)548 pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)549 pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)550 ENDDO551 ENDDO552 553 ENDIF554 555 !101 continue556 c-----------------------------------------------------------------------557 c On incremente les tendances physiques de la temperature du sol:558 c ---------------------------------------------------------------559 ! WRITE(55,*)"tsurf",maxval(tsurf(:)),minval(tsurf(:)),it560 561 DO j=jj_begin-offset,jj_end+offset562 DO i=ii_begin-offset,ii_end+offset563 ig=(j-1)*iim+i564 tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig)565 ENDDO566 ENDDO567 c-----------------------------------------------------------------------568 c soil temperatures:569 c --------------------570 571 IF (callsoil) THEN572 CALL soil(ngrid,nsoilmx,.false.,inertie,573 s ptimestep,tsurf,tsoil,capcal,fluxgrd)574 IF(lverbose) THEN575 PRINT*,'Surface Heat capacity,conduction Flux, Ts,576 s dTs, dt'577 PRINT*,capcal(igout),fluxgrd(igout),tsurf(igout),578 s zdtsrf(igout),ptimestep579 ENDIF580 ENDIF581 582 c-----------------------------------------------------------------------583 c sorties:584 c --------585 if(zday.GT.zday_last+period_sort) then586 zday_last=zday587 c Ecriture/extension de la coordonnee temps588 589 do ig=1,ngrid590 zpmer(ig)=pplev(ig,1)*exp(pphi(ig,1)/(kappa*cpp*285.))591 enddo592 endif593 c-----------------------------------------------------------------------594 IF(lastcall) THEN595 PRINT*,'Ecriture du fichier de reinitialiastion de la physique'596 ENDIF597 598 icount=icount+1599 ! RETURN600 END SUBROUTINE phyparam_lmd601 END MODULE PHY -
codes/icosagcm/trunk/src/physics.f90
r178 r186 2 2 3 3 CHARACTER(LEN=255) :: physics_type="automatic" 4 4 !$OMP THREADPRIVATE(physics_type) 5 5 6 6 CONTAINS … … 11 11 USE icosa 12 12 USE physics_dcmip_mod,init_physics_dcmip=>init_physics 13 USE physics_dry_mod13 ! USE physics_dry_mod 14 14 IMPLICIT NONE 15 15 … … 32 32 33 33 CASE ('dry') 34 CALL init_physics_dry34 ! CALL init_physics_dry 35 35 36 36 CASE DEFAULT … … 44 44 SUBROUTINE physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q) 45 45 USE icosa 46 USE physics_dry_mod46 ! USE physics_dry_mod 47 47 USE physics_dcmip_mod, physics_dcmip=>physics 48 48 USE etat0_mod … … 67 67 ! CALL transfert_request(f_ue,req_e1_vect) 68 68 CALL held_suarez(f_ps,f_theta_rhodz,f_ue) 69 CASE DEFAULT 70 ! PRINT*,"NO PHYSICAL PACAKAGE USED" ! FIXME MPI 69 71 END SELECT 70 72 … … 73 75 74 76 CASE ('dry') 75 CALL physics_dry(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)77 ! CALL physics_dry(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q) 76 78 77 79 CASE DEFAULT -
codes/icosagcm/trunk/src/physics_dcmip.f90
r146 r186 1 1 MODULE physics_dcmip_mod 2 2 USE ICOSA 3 3 PRIVATE 4 4 5 5 INTEGER,SAVE :: testcase 6 !$OMP THREADPRIVATE(testcase) 7 6 8 PUBLIC init_physics, physics 7 9 TYPE(t_field),POINTER :: f_out_i(:) … … 42 44 CALL transfert_request(f_ue,req_e1_vect) 43 45 DO ind=1,ndomain 46 IF (.NOT. assigned_domain(ind)) CYCLE 44 47 CALL swap_dimensions(ind) 45 48 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/physics_dry.f90
r149 r186 1 MODULE physics_dry_mod2 USE ICOSA3 PUBLIC init_physics_dry, physics_dry4 5 CONTAINS6 7 SUBROUTINE init_physics_dry8 USE ICOSA !sarvesh9 USE time_mod !sarvesh10 USE dimphys_mod11 USE RADIATION12 13 IMPLICIT NONE14 INTEGER::i,j,ij15 !-------------------------------------- ORBITAL PARAMETER----16 periheli=150.17 CALL getin('periheli', periheli)18 aphelie=150.19 CALL getin('aphelie',aphelie)20 coefir=0.0821 CALL getin('coefir',coefir)22 coefvis=0.9923 CALL getin('coefvis',coefvis)24 obliquit=0.025 CALL getin('obliquit',obliquit)26 peri_day=0.27 CALL getin('peri_day',peri_day)28 year_day=360.29 CALL getin('year_day',year_day)30 callrad=.true.31 CALL getin('callrad', callrad)32 calldifv=.true.33 CALL getin('calldifv', calldifv)34 calladj=.true.35 CALL getin('calladj', calladj)36 callcond=.true.37 callsoil=.true.38 CALL getin('callsoil',callsoil)39 season=.true.40 CALL getin('season',season)41 diurnal=.true.42 CALL getin('diurnal',diurnal)43 lverbose=.false.44 CALL getin('lverbose',lverbose)45 period_sort=1.46 CALL getin('period_sort',period_sort)47 ! ptimestep=dt48 ! CALL getin('ptimestep',ptimestep)49 50 print*,'Activation de la physique:'51 print*,' Rayonnement ',callrad52 print*,' Diffusion verticale turbulente ', calldifv53 print*,' Ajustement convectif ',calladj54 print*,' Sol ',callsoil55 print*,' Cycle diurne ',diurnal56 ! choice of the frequency of the computation of radiations57 IF(diurnal) THEN58 iradia=NINT(daysec/(20.*dt))59 ELSE60 iradia=NINT(daysec/(4.*dt))61 ENDIF62 iradia=163 64 ngridmx=iim*jjm ; nlayermx=llm65 offset=halo66 67 ALLOCATE(albedo(ngridmx));ALLOCATE(emissiv(ngridmx))68 ALLOCATE(inertie(ngridmx));ALLOCATE(z0(ngridmx))69 ALLOCATE(rnatur(ngridmx));ALLOCATE(tsurf(ngridmx))70 ALLOCATE(tsoil(ngridmx,nlayermx));ALLOCATE(fluxgrd(ngridmx))71 ALLOCATE(fluxrad(ngridmx));ALLOCATE(dtrad(ngridmx,llm+1))72 ALLOCATE(q2(ngridmx,llm+1));ALLOCATE(q2l(ngridmx,llm+1))73 ALLOCATE(capcal(ngridmx))74 75 CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)76 77 PRINT*,'unjours',daysec78 PRINT*,'The radiative transfer is computed each ', iradia,' physical time-step or each ', &79 iradia*dt,' seconds'80 END SUBROUTINE init_physics_dry81 82 83 SUBROUTINE physics_dry( it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)84 USE icosa85 IMPLICIT NONE86 INTEGER,INTENT(IN) :: it87 REAL(rstd),INTENT(IN) :: jD_cur,jH_cur88 TYPE(t_field),POINTER :: f_phis(:)89 TYPE(t_field),POINTER :: f_ps(:)90 TYPE(t_field),POINTER :: f_theta_rhodz(:)91 TYPE(t_field),POINTER :: f_ue(:)92 TYPE(t_field),POINTER :: f_q(:)93 94 REAL(rstd),POINTER :: phis(:)95 REAL(rstd),POINTER :: ps(:)96 REAL(rstd),POINTER :: theta_rhodz(:,:)97 REAL(rstd),POINTER :: ue(:,:)98 REAL(rstd),POINTER :: q(:,:,:)99 ! REAL(rstd),POINTER :: precl(:)100 INTEGER :: ind101 ! LOGICAL:: firstcall,lastcall102 103 CALL transfert_request(f_ue,req_e1_vect)104 CALL transfert_request(f_theta_rhodz,req_i1)105 106 DO ind=1,ndomain107 CALL swap_dimensions(ind)108 CALL swap_geometry(ind)109 phis=f_phis(ind)110 ps=f_ps(ind)111 theta_rhodz=f_theta_rhodz(ind)112 ue=f_ue(ind)113 q=f_q(ind)114 ! out_i=f_out_i(ind)115 ! precl=f_precl(ind)116 ! print*,"====================================ind",ind,"----------it",it117 CALL compute_physics_dry(it,jD_cur,jH_cur,phis, ps, theta_rhodz, ue, q(:,:,1))118 ENDDO119 120 ! CALL writefield("out_i",f_out_i)121 122 ! IF (mod(it,itau_out)==0 ) THEN123 ! CALL writefield("precl",f_precl)124 ! ENDIF125 126 END SUBROUTINE physics_dry127 128 SUBROUTINE compute_physics_dry(it,jD_cur,jH_cur,phis, ps, theta_rhodz, ue, q)129 USE icosa130 USE pression_mod131 USE exner_mod132 USE theta2theta_rhodz_mod133 USE geopotential_mod134 USE wind_mod135 USE PHY136 137 IMPLICIT NONE138 INTEGER::it139 REAL(rstd) :: jD_cur140 REAL(rstd) :: jH_cur141 REAL(rstd) :: phis(iim*jjm)142 REAL(rstd) :: ps(iim*jjm)143 REAL(rstd) :: theta_rhodz(iim*jjm,llm)144 REAL(rstd) :: ue(3*iim*jjm,llm)145 REAL(rstd) :: q(iim*jjm,llm)146 ! REAL(rstd) :: precl(iim*jjm)147 148 REAL(rstd) :: p(iim*jjm,llm+1)149 REAL(rstd) :: pks(iim*jjm)150 REAL(rstd) :: pk(iim*jjm,llm)151 REAL(rstd) :: phi(iim*jjm,llm)152 REAL(rstd) :: T(iim*jjm,llm)153 REAL(rstd) :: Tfi(iim*jjm,llm)154 REAL(rstd) :: theta(iim*jjm,llm)155 156 REAL(rstd) :: uc(iim*jjm,3,llm)157 REAL(rstd) :: u(iim*jjm,llm)158 REAL(rstd) :: v(iim*jjm,llm)159 REAL(rstd) :: ufi(iim*jjm,llm)160 REAL(rstd) :: vfi(iim*jjm,llm)161 REAL(rstd) :: qfi(iim*jjm,llm)162 REAL(rstd) :: utemp(iim*jjm,llm)163 REAL(rstd) :: vtemp(iim*jjm,llm)164 REAL(rstd) :: lat(iim*jjm)165 REAL(rstd) :: lon(iim*jjm)166 REAL(rstd) :: pmid(iim*jjm,llm)167 REAL(rstd) :: pint(iim*jjm,llm+1)168 REAL(rstd) :: pdel(iim*jjm,llm)169 REAL(rstd) :: plev(iim*jjm,llm+1),play(iim*jjm,llm)170 REAL(rstd) :: pkbycp171 INTEGER :: i,j,l,ij,ig172 173 !-------------------174 ! LOGICAL:: firstcall,lastcall175 REAL(rstd) :: dufi(iim*jjm,llm)176 REAL(rstd) :: dvfi(iim*jjm,llm)177 REAL(rstd) :: dTfi(iim*jjm,llm)178 REAL(rstd) :: dpsfi(iim*jjm)179 REAL(rstd) :: dqfi(iim*jjm,llm)180 ! PRINT *,'Entering in LMD SIMPLE physics'181 182 183 offset=halo184 CALL compute_pression(ps,p,halo)185 CALL compute_exner(ps,p,pks,pk,halo)186 CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,halo)187 CALL compute_geopotential(phis,pks,pk,theta,phi,halo)188 CALL compute_theta_rhodz2temperature(ps,theta_rhodz,T,halo)189 CALL compute_wind_centered(ue,uc)190 CALL compute_wind_centered_lonlat_compound(uc, u, v)191 192 DO j=jj_begin-offset,jj_end+offset193 DO i=ii_begin-offset,ii_end+offset194 ij=(j-1)*iim+i195 CALL xyz2lonlat(xyz_i(ij,:),lon(ij),lat(ij))196 ENDDO197 ENDDO198 199 DO l=1,llm200 DO j=jj_begin-offset,jj_end+offset201 DO i=ii_begin-offset,ii_end+offset202 ij=(j-1)*iim+i203 ! Tfi(ij,l)=T(ij,l)204 ! ufi(ij,l)=u(ij,l)205 ! vfi(ij,l)=v(ij,l)206 ! qfi(ij,l)=q(ij,l)207 dTfi(ij,l)=0.0208 dufi(ij,l)=0.0209 dvfi(ij,l)=0.0210 dqfi(ij,l)=0.0211 ENDDO212 ENDDO213 ENDDO214 plev(:,:) = p(:,:)215 dpsfi=0.0216 217 DO l=1,llm218 DO j=jj_begin-offset,jj_end+offset219 DO i=ii_begin-offset,ii_end+offset220 ij=(j-1)*iim+i221 pkbycp=pk(ij,l)/cpp222 play(ij,l)=preff*pkbycp**(1./kappa)223 ENDDO224 ENDDO225 ENDDO226 227 228 CALL phyparam_lmd(it,iim*jjm,llm,1,dt,lat,lon,jD_cur,jH_cur, &229 plev,play,phi,phis,u,v,T,q,dufi,dvfi,dTfi,dqfi,dpsfi)230 231 CALL ADDFI(u,v,T,q,ps,dufi,dvfi,dTfi,dqfi,dpsfi)232 233 ! CALL SARCHECKF(llm)234 ! print*,"plev",(maxval(plev(:,l)),l=1,llm+1)235 236 ! CALL phyparam_lmd(it,iim*jjm,llm,1,dt,lat,lon,jD_cur,jH_cur, &237 ! plev,play,phi,phis,ufi,vfi,Tfi,qfi,dufi,dvfi,dTfi,dqfi,dpsfi)238 239 ! Print*,"going ADD FI",it240 ! CALL ADDFI(ufi,vfi,Tfi,qfi,ps,dufi,dvfi,dTfi,dqfi,dpsfi)241 242 ! WRITE(11,*)"ducovfi",maxval(dufi),minval(dufi),it243 ! WRITE(11,*)"ucovfi",maxval(ufi),minval(ufi)244 ! WRITE(11,*)"dtetafi",maxval(dTfi),minval(dTfi)245 246 !=============================================247 ! go to 1234248 DO l=1,llm249 DO j=jj_begin-offset,jj_end+offset250 DO i=ii_begin-offset,ii_end+offset251 ij=(j-1)*iim+i252 uc(ij,:,l)=(dufi(ij,l)*elon_i(ij,:)+dvfi(ij,l)*elat_i(ij,:))*dt253 ENDDO254 ENDDO255 ENDDO256 257 258 DO l=1,llm259 DO j=jj_begin-offset,jj_end+offset260 DO i=ii_begin-offset,ii_end+offset261 ij=(j-1)*iim+i262 ue(ij+u_right,l)=ue(ij+u_right,l)+sum( 0.5*(uc(ij,:,l) + uc(ij+t_right,:,l))*ep_e(ij+u_right,:) )263 ue(ij+u_lup,l)=ue(ij+u_lup,l)+sum( 0.5*(uc(ij,:,l) + uc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )264 ue(ij+u_ldown,l)=ue(ij+u_ldown,l)+sum( 0.5*(uc(ij,:,l) + uc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )265 ENDDO266 ENDDO267 ENDDO268 !1234 continue269 270 ! CALL compute_temperature2theta_rhodz(ps,Tfi,theta_rhodz,halo)271 CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,halo)272 273 ! WRITE(13,*)"tetafi",maxval(theta_rhodz),minval(theta_rhodz)274 RETURN275 END SUBROUTINE compute_physics_dry276 277 SUBROUTINE addfi(ufi,vfi,Tfi,qfi,ps,dufi,dvfi,dTfi,dqfi,dpsfi)278 USE ICOSA279 IMPLICIT NONE280 REAL(rstd) :: dufi(iim*jjm,llm)281 REAL(rstd) :: dvfi(iim*jjm,llm)282 REAL(rstd) :: dTfi(iim*jjm,llm)283 REAL(rstd) :: dpsfi(iim*jjm)284 REAL(rstd) :: dqfi(iim*jjm,llm)285 REAL(rstd) :: ufi(iim*jjm,llm)286 REAL(rstd) :: vfi(iim*jjm,llm)287 REAL(rstd) :: qfi(iim*jjm,llm)288 REAL(rstd) :: ps(iim*jjm)289 REAL(rstd) :: Tfi(iim*jjm,llm)290 INTEGER::i,j,l,ij,offset291 offset=halo292 293 DO l=1,llm294 DO j=jj_begin-offset,jj_end+offset295 DO i=ii_begin-offset,ii_end+offset296 ij=(j-1)*iim+i297 Tfi(ij,l)=Tfi(ij,l)+dTfi(ij,l)*dt298 ufi(ij,l)=ufi(ij,l)+dufi(ij,l)*dt299 vfi(ij,l)=vfi(ij,l)+dvfi(ij,l)*dt300 qfi(ij,l)=qfi(ij,l)+dqfi(ij,l)*dt301 END DO302 END DO303 END DO304 ps(:)=ps(:) + dpsfi(:)*dt305 END SUBROUTINE addfi306 307 END MODULE physics_dry_mod -
codes/icosagcm/trunk/src/pression.f90
r166 r186 14 14 15 15 DO ind=1,ndomain 16 IF (.NOT. assigned_domain(ind)) CYCLE 16 17 CALL swap_dimensions(ind) 17 18 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/radiation_mod.F
r149 r186 1 MODULE RADIATION2 USE ICOSA3 USE dimphys_mod4 ! USE PHY5 contains6 7 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%8 SUBROUTINE zerophys(n,x)9 IMPLICIT NONE10 INTEGER n11 REAL x(n)12 INTEGER i13 14 DO i=1,n15 x(i)=0.16 ENDDO17 RETURN18 END subroutine zerophys19 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%20 21 SUBROUTINE solarlong(pday,psollong)22 IMPLICIT NONE23 c=======================================================================24 c25 c Objet:26 c ------27 c28 c Calcul de la distance soleil-planete et de la declinaison29 c en fonction du jour de l'annee.30 c31 c32 c Methode:33 c --------34 c35 c Calcul complet de l'elipse36 c37 c Interface:38 c ----------39 c40 c Uncommon comprenant les parametres orbitaux.41 c42 c Arguments:43 c ----------44 c45 c Input:46 c ------47 c pday jour de l'annee (le jour 0 correspondant a l'equinoxe)48 c lwrite clef logique pour sorties de controle49 c50 c Output:51 c -------52 c pdist_sol distance entre le soleil et la planete53 c ( en unite astronomique pour utiliser la constante54 c solaire terrestre 1370 Wm-2 )55 c pdecli declinaison ( en radians )56 c57 c=======================================================================58 c-----------------------------------------------------------------------59 c Declarations:60 c -------------61 62 c arguments:63 c ----------64 65 REAL pday,pdist_sol,pdecli,psollong66 LOGICAL lwrite67 68 c Local:69 c ------70 71 REAL zanom,xref,zx0,zdx,zteta,zz72 INTEGER iter73 74 75 c-----------------------------------------------------------------------76 c calcul de l'angle polaire et de la distance au soleil :77 c -------------------------------------------------------78 79 c calcul de l'zanomalie moyenne80 81 zz=(pday-peri_day)/year_day82 zanom=2.*pi*(zz-nint(zz))83 xref=abs(zanom)84 85 c resolution de lequation horaire zx0 - e * sin (zx0) = xref86 c methode de Newton87 88 zx0=xref+e_elips*sin(xref)89 DO 110 iter=1,1090 zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))91 if(abs(zdx).le.(1.e-7)) goto 12092 zx0=zx0+zdx93 110 continue94 120 continue95 zx0=zx0+zdx96 if(zanom.lt.0.) zx0=-zx097 98 c zteta est la longitude solaire99 100 zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))101 102 psollong=zteta-timeperi103 104 IF(psollong.LT.0.) psollong=psollong+2.*pi105 IF(psollong.GT.2.*pi) psollong=psollong-2.*pi106 107 RETURN108 END SUBROUTINE solarlong109 110 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%111 112 SUBROUTINE orbite(pls,pdist_sol,pdecli)113 IMPLICIT NONE114 c====================================================================115 c116 c Objet:117 c ------118 c119 c Distance from sun and declimation as a function of the solar120 c longitude Ls121 c122 c Interface:123 c ----------124 c Arguments:125 c ----------126 c127 c Input:128 c ------129 c pls Ls130 c131 c Output:132 c -------133 c pdist_sol Distance Sun-Planet in UA134 c pdecli declinaison ( en radians )135 c136 c====================================================================137 c-----------------------------------------------------------------------138 c Declarations:139 c -------------140 c arguments:141 c ----------142 143 REAL pday,pdist_sol,pdecli,pls144 145 c-----------------------------------------------------------------------146 147 c Distance Sun-Planet148 149 pdist_sol=p_elips/(1.+e_elips*cos(pls+timeperi))150 151 c Solar declination152 153 pdecli= asin (sin(pls)*sin(obliquit*pi/180.))154 155 c-----------------------------------------------------------------------156 c sorties eventuelles:157 c ---------------------158 159 END SUBROUTINE orbite160 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%161 162 SUBROUTINE iniorbit163 $ (paphelie,pperiheli,pyear_day,pperi_day,pobliq)164 IMPLICIT NONE165 c=====================================================================166 c167 c Auteur:168 c -------169 c Frederic Hourdin 22 Fevrier 1991170 c171 c Objet:172 c ------173 c Initialisation du sous programme orbite qui calcule174 c a une date donnee de l'annee de duree year_day commencant175 c a l'equinoxe de printemps et dont le perihelie se situe176 c a la date peri_day, la distance au soleil et la declinaison.177 c178 c Interface:179 c ----------180 c - Doit etre appele avant d'utiliser orbite.181 c - initialise le common comorbit182 c183 c Arguments:184 c ----------185 c186 c Input:187 c ------188 c aphelie \ aphelie et perihelie de l'orbite189 c periheli / en millions de kilometres.190 c191 c=====================================================================192 c Declarations:193 c -------------194 c Arguments:195 c ----------196 197 REAL paphelie,pperiheli,pyear_day,pperi_day,pobliq198 199 c Local:200 c ------201 202 REAL zxref,zanom,zz,zx0,zdx203 INTEGER iter204 205 c'-----------------------------------------------------------------------206 207 aphelie =paphelie208 periheli=pperiheli209 year_day=pyear_day210 obliquit=pobliq211 peri_day=pperi_day212 213 PRINT*,'Perihelie en Mkm ',periheli214 PRINT*,'Aphelise en Mkm ',aphelie215 PRINT*,'obliquite en degres :',obliquit216 unitastr=149.597927217 e_elips=(aphelie-periheli)/(periheli+aphelie)218 p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr219 220 print*,'e_elips',e_elips221 print*,'p_elips',p_elips222 223 c-----------------------------------------------------------------------224 c calcul de l'angle polaire et de la distance au soleil :225 c -------------------------------------------------------226 227 c calcul de l'zanomalie moyenne228 229 zz=(year_day-pperi_day)/year_day230 zanom=2.*pi*(zz-nint(zz))231 zxref=abs(zanom)232 PRINT*,'zanom ',zanom233 234 c resolution de lequation horaire zx0 - e * sin (zx0) = zxref235 c methode de Newton236 237 zx0=zxref+e_elips*sin(zxref)238 DO 110 iter=1,100239 zdx=-(zx0-e_elips*sin(zx0)-zxref)/(1.-e_elips*cos(zx0))240 if(abs(zdx).le.(1.e-12)) goto 120241 zx0=zx0+zdx242 110 continue243 120 continue244 zx0=zx0+zdx245 if(zanom.lt.0.) zx0=-zx0246 PRINT*,'zx0 ',zx0247 248 c zteta est la longitude solaire249 250 timeperi=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))251 PRINT*,'longitude solaire du perihelie timeperi = ',timeperi252 253 RETURN254 END SUBROUTINE iniorbit255 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%256 257 SUBROUTINE mucorr(npts,pdeclin, plat, pmu, pfract,phaut,prad)258 IMPLICIT NONE259 260 c=======================================================================261 c262 c Calcul of equivalent solar angle and and fraction of day whithout263 c diurnal cycle.264 c265 c parmeters :266 c -----------267 c268 c Input :269 c -------270 c npts number of points271 c pdeclin solar declinaison272 c plat(npts) latitude273 c phaut hauteur typique de l'atmosphere274 c prad rayon planetaire '275 c276 c Output :277 c --------278 c pmu(npts) equivalent cosinus of the solar angle279 c pfract(npts) fractionnal day280 c281 c=======================================================================282 283 c-----------------------------------------------------------------------284 c285 c 0. Declarations :286 c -----------------287 288 c Arguments :289 c -----------290 INTEGER npts291 REAL plat(npts),pmu(npts), pfract(npts)292 REAL phaut,prad,pdeclin293 c294 c Local variables :295 c -----------------296 INTEGER j,i,ij,ig297 REAL pi298 REAL z,cz,sz,tz,phi,cphi,sphi,tphi299 REAL ap,a,t,b300 REAL alph301 302 c-----------------------------------------------------------------------303 304 !print*,'npts,pdeclin'305 !print*,npts,pdeclin306 pi = 4. * atan(1.0)307 !print*,'PI=',pi308 pi=2.*asin(1.)309 z = pdeclin310 cz = cos (z)311 sz = sin (z)312 !print*,'cz,sz',cz,sz313 314 ! DO j=jj_begin-offset,jj_end+offset315 ! DO i=ii_begin-offset,ii_end+offset316 ! ig=(j-1)*iim+i317 318 DO ig=1,npts319 phi = plat(ig)320 cphi = cos(phi)321 if (cphi.le.1.e-9) cphi=1.e-9322 sphi = sin(phi)323 tphi = sphi / cphi324 b = cphi * cz325 t = -tphi * sz / cz326 a = 1.0 - t*t327 ap = a328 IF(t.eq.0.) then329 t=0.5*pi330 ELSE331 IF (a.lt.0.) a = 0.332 t = sqrt(a) / t333 IF (t.lt.0.) then334 t = -atan (-t) + pi335 ELSE336 t = atan(t)337 ENDIF338 ENDIF339 340 pmu(ig) = (sphi*sz*t) / pi + b*sin(t)/pi341 pfract(ig) = t / pi342 IF (ap .lt.0.) then343 pmu(ig) = sphi * sz344 pfract(ig) = 1.0345 ENDIF346 IF (pmu(ig).le.0.0) pmu(ig) = 0.0347 pmu(ig) = pmu(ig) / pfract(ig)348 IF (pmu(ig).eq.0.) pfract(ig) = 0.349 ENDDO350 ! ENDDO351 c-----------------------------------------------------------------------352 c correction de rotondite:353 c ------------------------354 355 alph=phaut/prad356 ! DO j=jj_begin-offset,jj_end+offset357 ! DO i=ii_begin-offset,ii_end+offset358 ! ig=(j-1)*iim+i359 DO ig = 1,npts360 pmu(ig)=sqrt(1224.*pmu(ig)*pmu(ig)+1.)/35.361 362 c $ (sqrt(alph*alph*pmu(ig)*pmu(ig)+2.*alph+1.)-alph*pmu(ig))363 ENDDO364 ! ENDDO365 366 RETURN367 END SUBROUTINE mucorr368 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%369 370 SUBROUTINE sw(ngrid,nlayer,ldiurn,371 $ coefvis,albedo,372 $ plevel,ps_rad,pmu,pfract,psolarf0,373 $ fsrfvis,dtsw,374 $ lwrite)375 IMPLICIT NONE376 c=======================================================================377 c378 c Rayonnement solaire en atmosphere non diffusante avec un379 c coefficient d'absoprption gris.380 c'381 c=======================================================================382 c383 c declarations:384 c -------------385 c arguments:386 c ----------387 c388 INTEGER ngrid,nlayer,i,j,ij389 REAL albedo(ngrid),coefvis390 REAL pmu(ngrid),pfract(ngrid)391 REAL plevel(ngrid,nlayer+1),ps_rad392 REAL psolarf0393 REAL fsrfvis(ngrid),dtsw(ngrid,nlayer)394 LOGICAL lwrite,ldiurn395 c396 c variables locales:397 c ------------------398 c399 400 REAL zalb(ngrid),zmu(ngrid),zfract(ngrid)401 REAL zplev(ngrid,nlayer+1)402 REAL zflux(ngrid),zdtsw(ngrid,nlayer)403 404 INTEGER ig,l,nlevel,ncount,igout405 INTEGER,DIMENSION(ngrid)::index406 REAL ztrdir(ngrid,nlayer+1),ztrref(ngrid,nlayer+1)407 REAL zfsrfref(ngrid)408 REAL z1(ngrid)409 REAL zu(ngrid,nlayer+1)410 REAL tau0411 412 EXTERNAL SSUM413 EXTERNAL ismax414 REAL ismax415 416 LOGICAL firstcall417 SAVE firstcall418 DATA firstcall/.true./419 420 c-----------------------------------------------------------------------421 c 1. initialisations:422 c -------------------423 424 425 426 IF (firstcall) THEN427 428 IF (ngrid.NE.ngrid) THEN429 PRINT*,'STOP in inifis'430 PRINT*,'Probleme de dimenesions :'431 PRINT*,'ngrid = ',ngrid432 PRINT*,'ngrid = ',ngrid433 STOP434 ENDIF435 436 437 IF (nlayer.NE.llm) THEN438 PRINT*,'STOP in inifis'439 PRINT*,'Probleme de dimenesions :'440 PRINT*,'nlayer = ',nlayer441 PRINT*,'llm = ',llm442 STOP443 ENDIF444 445 ENDIF446 447 nlevel=nlayer+1448 c-----------------------------------------------------------------------449 c Definitions des tableaux locaux pour les points ensoleilles:450 c ------------------------------------------------------------451 452 IF (ldiurn) THEN453 ncount=0454 DO j=jj_begin-offset,jj_end+offset455 DO i=ii_begin-offset,ii_end+offset456 ig=(j-1)*iim+i457 index(ig)=0458 ENDDO459 ENDDO460 DO j=jj_begin-offset,jj_end+offset461 DO i=ii_begin-offset,ii_end+offset462 ig=(j-1)*iim+i463 IF(pfract(ig).GT.1.e-6) THEN464 ncount=ncount+1465 index(ncount)=ig466 ENDIF467 ENDDO468 ENDDO469 ! SARVESH CHANGED NCOUNT TO NGRID TO BE CONSISTENT ???470 471 CALL monGATHER(ncount,zfract,pfract,index)472 CALL monGATHER(ncount,zmu,pmu,index)473 CALL monGATHER(ncount,zalb,albedo,index)474 DO l=1,nlevel475 CALL monGATHER(ncount,zplev(1,l),plevel(1,l),index)476 ENDDO477 ELSE478 ncount=ngrid479 zfract(:)=pfract(:)480 zmu(:)=pmu(:)481 zalb(:)=albedo(:)482 zplev(:,:)=plevel(:,:)483 ENDIF484 485 c-----------------------------------------------------------------------486 c calcul des profondeurs optiques integres depuis p=0:487 c ----------------------------------------------------488 489 tau0=-.5*log(coefvis)490 491 c calcul de la partie homogene de l'opacite492 c'493 tau0=tau0/ps_rad494 495 496 DO l=1,nlayer+1497 DO j=jj_begin-offset,jj_end+offset498 DO i=ii_begin-offset,ii_end+offset499 ig=(j-1)*iim+i500 zu(ig,l)=tau0*zplev(ig,l)501 ENDDO502 ENDDO503 ENDDO504 505 c-----------------------------------------------------------------------506 c 2. calcul de la transmission depuis le sommet de l'atmosphere:507 c' -----------------------------------------------------------508 509 DO l=1,nlevel510 DO j=jj_begin-offset,jj_end+offset511 DO i=ii_begin-offset,ii_end+offset512 ig=(j-1)*iim+i513 ztrdir(ig,l)=exp(-zu(ig,l)/zmu(ig))514 ENDDO515 ENDDO516 ENDDO517 518 IF (lwrite) THEN519 igout=ncount/2+1520 PRINT*521 PRINT*,'Diagnostique des transmission dans le spectre solaire'522 PRINT*,'zfract, zmu, zalb'523 PRINT*,zfract(igout), zmu(igout), zalb(igout)524 PRINT*,'Pression, quantite d abs, transmission'525 DO l=1,nlayer+1526 PRINT*,zplev(igout,l),zu(igout,l),ztrdir(igout,l)527 ENDDO528 ENDIF529 530 c-----------------------------------------------------------------------531 c 3. taux de chauffage, ray. solaire direct:532 c ------------------------------------------533 534 DO l=1,nlayer535 DO j=jj_begin-offset,jj_end+offset536 DO i=ii_begin-offset,ii_end+offset537 ig=(j-1)*iim+i538 zdtsw(ig,l)=g*psolarf0*zfract(ig)*zmu(ig)*539 $ (ztrdir(ig,l+1)-ztrdir(ig,l))/540 $ (cpp*(zplev(ig,l)-zplev(ig,l+1)))541 ENDDO542 ENDDO543 ENDDO544 IF (lwrite) THEN545 PRINT*546 PRINT*,'Diagnostique des taux de chauffage solaires:'547 PRINT*,' 1 taux de chauffage lie au ray. solaire direct'548 DO l=1,nlayer549 PRINT*,zdtsw(igout,l)550 ENDDO551 ENDIF552 553 554 c-----------------------------------------------------------------------555 c 4. calcul du flux solaire arrivant sur le sol:556 c ----------------------------------------------557 558 DO j=jj_begin-offset,jj_end+offset559 DO i=ii_begin-offset,ii_end+offset560 ig=(j-1)*iim+i561 z1(ig)=zfract(ig)*zmu(ig)*psolarf0*ztrdir(ig,1)562 zflux(ig)=(1.-zalb(ig))*z1(ig)563 zfsrfref(ig)= zalb(ig)*z1(ig)564 ENDDO565 ENDDO566 IF (lwrite) THEN567 PRINT*568 PRINT*,'Diagnostique des taux de chauffage solaires:'569 PRINT*,' 2 flux solaire net incident sur le sol'570 PRINT*,zflux(igout)571 ENDIF572 573 574 c-----------------------------------------------------------------------575 c 5.calcul des traansmissions depuis le sol, cas diffus:576 c ------------------------------------------------------577 578 DO l=1,nlevel579 DO j=jj_begin-offset,jj_end+offset580 DO i=ii_begin-offset,ii_end+offset581 ig=(j-1)*iim+i582 ztrref(ig,l)=exp(-(zu(ig,1)-zu(ig,l))*1.66)583 ENDDO584 ENDDO585 ENDDO586 587 IF (lwrite) THEN588 PRINT*589 PRINT*,'Diagnostique des taux de chauffage solaires'590 PRINT*,' 3 transmission avec les sol'591 PRINT*,'niveau transmission'592 DO l=1,nlevel593 PRINT*,l,ztrref(igout,l)594 ENDDO595 ENDIF596 597 c-----------------------------------------------------------------------598 c 6.ajout a l'echauffement de la contribution du ray. sol. reflechit:599 c' -------------------------------------------------------------------600 601 DO l=1,nlayer602 DO j=jj_begin-offset,jj_end+offset603 DO i=ii_begin-offset,ii_end+offset604 ig=(j-1)*iim+i605 zdtsw(ig,l)=zdtsw(ig,l)+606 $ g*zfsrfref(ig)*(ztrref(ig,l+1)-ztrref(ig,l))/607 $ (cpp*(zplev(ig,l+1)-zplev(ig,l)))608 ENDDO609 ENDDO610 ENDDO611 612 c-----------------------------------------------------------------------613 c 10. sorites eventuelles:614 c ------------------------615 616 IF (lwrite) THEN617 PRINT*618 PRINT*,'Diagnostique des taux de chauffage solaires:'619 PRINT*,' 3 taux de chauffage total'620 DO l=1,nlayer621 PRINT*,zdtsw(igout,l)622 ENDDO623 ENDIF624 625 IF (ldiurn) THEN626 CALL zerophys(ngrid,fsrfvis)627 CALL monscatter(ncount,fsrfvis,index,zflux)628 CALL zerophys(ngrid*nlayer,dtsw)629 DO l=1,nlayer630 CALL monscatter(ncount,dtsw(1,l),index,zdtsw(1,l))631 ENDDO632 ELSE633 !print*,'NOT DIURNE'634 fsrfvis(:)=zflux(:)635 dtsw(:,:)=zdtsw(:,:)636 ENDIF637 638 RETURN639 END SUBROUTINE sw640 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%641 SUBROUTINE lw(ngrid,nlayer,coefir,emissiv,642 $ pp,ps_rad,ptsurf,pt,643 $ pfluxir,pdtlw,644 $ lwrite)645 646 IMPLICIT NONE647 c=======================================================================648 c649 c calcul de l'evolution de la temperature sous l'effet du rayonnement650 c infra-rouge.651 c Pour simplifier, les transmissions sont precalculees et ne652 c dependent que de l'altitude.653 c654 c arguments:655 c ----------656 c'657 c entree:658 c -------659 c ngrid nombres de points de la grille horizontale660 c nlayer nombre de couches661 c ptsurf(ngrid) temperature de la surface662 c pt(ngrid,nlayer) temperature des couches663 c pp(ngrid,nlayer+1) pression entre les couches664 c lwrite variable logique pour sorties665 c666 c sortie:667 c -------668 c pdtlw(ngrid,nlayer) taux de refroidissement669 c pfluxir(ngrid) flux infrarouge sur le sol670 c671 c=======================================================================672 673 !c declarations:674 !c -------------675 !c arguments:676 !c' ----------677 678 INTEGER ngrid,nlayer679 REAL coefir,emissiv(ngrid),ps_rad680 REAL ptsurf(ngrid),pt(ngrid,nlayer),pp(ngrid,nlayer+1)681 REAL pdtlw(ngrid,nlayer),pfluxir(ngrid)682 LOGICAL lwrite683 684 c variables locales:685 c ------------------686 687 INTEGER nlevel,ilev,ig,i,il688 REAL zplanck(ngridmx,nlayermx+1),zcoef689 REAL zfluxup(ngridmx,nlayermx+1),zfluxdn(ngridmx,nlayermx+1)690 REAL zflux(ngridmx,nlayermx+1)691 REAL zlwtr1(ngridmx),zlwtr2(ngridmx)692 REAL zup(ngridmx,nlayermx+1),zdup(ngridmx)693 REAL stephan694 695 LOGICAL lstrong696 SAVE lstrong,stephan697 DATA lstrong/.true./698 c-----------------------------------------------------------------------699 c initialisations:700 c ----------------701 702 nlevel=nlayer+1703 stephan=5.67e-08704 705 706 pfluxir=0.0707 pdtlw=0.0708 !print*,"ngr,nlay,coefi",ngrid,nlayer,coefir709 c-----------------------------------------------------------------------710 c 2. calcul des quantites d'absorbants:711 c' -------------------------------------712 713 c absorption forte714 IF(lstrong) THEN715 DO ilev=1,nlevel716 DO ig=1,ngrid717 zup(ig,ilev)=pp(ig,ilev)*pp(ig,ilev)/(2.*g)718 ENDDO719 ENDDO720 IF(lwrite) THEN721 DO ilev=1,nlayer722 PRINT*,' up(',ilev,') = ',zup(ngrid/2+1,ilev)723 ENDDO724 ENDIF725 zcoef=-log(coefir)/sqrt(ps_rad*ps_rad/(2.*g))726 727 c absorption faible728 ELSE729 DO ilev=1,nlevel730 DO ig=1,ngrid731 zup(ig,ilev)=pp(ig,ilev)732 ENDDO733 ENDDO734 zcoef=-log(coefir)/ps_rad735 ENDIF736 737 738 c-----------------------------------------------------------------------739 c 2. calcul de la fonction de corps noir:740 c ---------------------------------------741 742 DO ilev=1,nlayer743 DO ig=1,ngrid744 zplanck(ig,ilev)=pt(ig,ilev)*pt(ig,ilev)745 zplanck(ig,ilev)=stephan*746 $ zplanck(ig,ilev)*zplanck(ig,ilev)747 ENDDO748 ENDDO749 750 c-----------------------------------------------------------------------751 c 4. flux descendants:752 c --------------------753 754 DO ilev=1,nlayer755 DO ig=1,ngrid756 zfluxdn(ig,ilev)=0.757 ENDDO758 DO ig=1,ngrid759 zdup(ig)=zup(ig,ilev)-zup(ig,nlevel)760 ENDDO761 CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)762 763 DO il=nlayer,ilev,-1764 zlwtr2(:)=zlwtr1(:)765 DO ig=1,ngrid766 zdup(ig)=zup(ig,ilev)-zup(ig,il)767 ENDDO768 CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)769 DO ig=1,ngrid770 zfluxdn(ig,ilev)=zfluxdn(ig,ilev)+771 $ zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))772 ENDDO773 ENDDO774 ENDDO775 776 DO ig=1,ngrid777 zfluxdn(ig,nlevel)=0.778 pfluxir(ig)=emissiv(ig)*zfluxdn(ig,1)779 ENDDO780 781 DO ig=1,ngrid782 zfluxup(ig,1)=ptsurf(ig)*ptsurf(ig)783 zfluxup(ig,1)=emissiv(ig)*stephan*zfluxup(ig,1)*zfluxup(ig,1)784 $ +(1.-emissiv(ig))*zfluxdn(ig,1)785 ENDDO786 787 c-----------------------------------------------------------------------788 c 3. flux montants:789 c ------------------790 791 DO ilev=1,nlayer792 DO ig=1,ngrid793 zdup(ig)=zup(ig,1)-zup(ig,ilev+1)794 ENDDO795 CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)796 DO ig=1,ngrid797 zfluxup(ig,ilev+1)=zfluxup(ig,1)*zlwtr1(ig)798 ENDDO799 DO il=1,ilev800 zlwtr2(:)=zlwtr1(:)801 DO ig=1,ngrid802 zdup(ig)=zup(ig,il+1)-zup(ig,ilev+1)803 ENDDO804 CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)805 DO ig=1,ngrid806 zfluxup(ig,ilev+1)=zfluxup(ig,ilev+1)+807 $ zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))808 ENDDO809 ENDDO810 811 ENDDO812 813 c-----------------------------------------------------------------------814 c 5. calcul des flux nets:815 c ------------------------816 817 DO ilev=1,nlevel818 DO ig=1,ngrid819 zflux(ig,ilev)=zfluxup(ig,ilev)-zfluxdn(ig,ilev)820 ENDDO821 ENDDO822 823 c-----------------------------------------------------------------------824 c 6. Calcul des taux de refroidissement:825 c --------------------------------------826 827 DO ilev=1,nlayer828 DO ig=1,ngrid829 pdtlw(ig,ilev)=(zflux(ig,ilev+1)-zflux(ig,ilev))*830 $ g/(cpp*(pp(ig,ilev+1)-pp(ig,ilev)))831 ENDDO832 ENDDO833 834 c-----------------------------------------------------------------------835 c 10. sorties eventuelles:836 c ------------------------837 838 IF (lwrite) THEN839 840 PRINT*841 PRINT*,'Diagnostique rayonnement thermique'842 PRINT*843 PRINT*,'temperature ',844 $ 'flux montant flux desc. taux de refroid.'845 i=ngrid/2+1846 WRITE(6,9000) ptsurf(i)847 DO ilev=1,nlayer848 WRITE(6,'(i4,4e18.4)') ilev,pt(i,ilev),849 $ zfluxup(i,ilev),zfluxdn(i,ilev),pdtlw(i,ilev)850 ENDDO851 WRITE(6,9000) zfluxup(i,nlevel),zfluxdn(i,nlevel)852 853 ENDIF854 855 c-----------------------------------------------------------------------856 857 RETURN858 9000 FORMAT(4e18.4)859 END SUBROUTINE lw860 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%861 862 subroutine solang ( kgrid,psilon,pcolon,psilat,pcolat,863 & ptim1,ptim2,ptim3,pmu0,pfract )864 IMPLICIT NONE865 866 C867 C**** *LW* - ORGANIZES THE LONGWAVE CALCULATIONS868 C869 C PURPOSE.870 C --------871 C CALCULATES THE SOLAR ANGLE FOR ALL THE POINTS OF THE GRID872 C ==== INPUTS ===873 C874 C PSILON(KGRID) : SINUS OF THE LONGITUDE875 C PCOLON(KGRID) : COSINUS OF THE LONGITUDE876 C PSILAT(KGRID) : SINUS OF THE LATITUDE877 C PCOLAT(KGRID) : COSINUS OF THE LATITUDE878 C PTIM1 : SIN(DECLI)879 C PTIM2 : COS(DECLI)*COS(TIME)880 C PTIM3 : SIN(DECLI)*SIN(TIME)881 C882 C ==== OUTPUTS ===883 C884 C PMU0 (KGRID) : SOLAR ANGLE885 C PFRACT(KGRID) : DAY FRACTION OF THE TIME INTERVAL886 C887 C IMPLICIT ARGUMENTS : NONE888 C --------------------889 C890 C METHOD.891 C -------892 C893 C EXTERNALS.894 C ----------895 C896 C NONE897 C898 C REFERENCE.899 C ----------900 C901 C RADIATIVE PROCESSES IN METEOROLOGIE AND CLIMATOLOGIE902 C PALTRIDGE AND PLATT903 C904 C AUTHOR.905 C -------906 C FREDERIC HOURDIN907 C908 C MODIFICATIONS.909 C --------------910 C ORIGINAL :90-01-14911 C 92-02-14 CALCULATIONS DONE THE ENTIER GRID (J.Polcher)912 C-----------------------------------------------------------------------913 C914 C ------------------------------------------------------------------915 916 C-----------------------------------------------------------------------917 C918 C* 0.1 ARGUMENTS919 C ---------920 C921 INTEGER kgrid922 REAL ptim1,ptim2,ptim3923 REAL psilon(kgrid),pcolon(kgrid),pmu0(kgrid),pfract(kgrid)924 REAL psilat(kgrid), pcolat(kgrid)925 C926 INTEGER jl,i,j927 REAL ztim1,ztim2,ztim3928 C------------------------------------------------------------------------929 C------------------------------------------------------------------------930 C---------------------------------------------------------------------931 C932 C--------------------------------------------------------------------933 C934 C* 1. INITIALISATION935 C --------------936 C937 !!!!!! 100 CONTINUE938 C939 DO j=jj_begin-offset,jj_end+offset940 DO i=ii_begin-offset,ii_end+offset941 jl=(j-1)*iim+i942 pmu0(jl)=0.943 pfract(jl)=0.944 ENDDO945 ENDDO946 !C947 !C* 1.1 COMPUTATION OF THE SOLAR ANGLE948 !C ------------------------------949 !C950 DO j=jj_begin-offset,jj_end+offset951 DO i=ii_begin-offset,ii_end+offset952 jl=(j-1)*iim+i953 ztim1=psilat(jl)*ptim1954 ztim2=pcolat(jl)*ptim2955 ztim3=pcolat(jl)*ptim3956 pmu0(jl)=ztim1+ztim2*pcolon(jl)+ztim3*psilon(jl)957 ENDDO958 ENDDO959 !C960 !C* 1.2 DISTINCTION BETWEEN DAY AND NIGHT961 !C ---------------------------------962 !C963 DO j=jj_begin-offset,jj_end+offset964 DO i=ii_begin-offset,ii_end+offset965 jl=(j-1)*iim+i966 IF (pmu0(jl).gt.0.) THEN967 pfract(jl)=1.968 ELSE969 pmu0(jl)=0.970 pfract(jl)=0.971 ENDIF972 ENDDO973 ENDDO974 RETURN975 END SUBROUTINE solang976 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%977 978 SUBROUTINE monGATHER(n,a,b,index)979 c980 IMPLICIT NONE981 C982 INTEGER n,ng,index(n),i,j,ij983 REAL a(n),b(n)984 c985 DO 100 i=1,n986 a(i)=b(index(i))987 100 CONTINUE988 989 ! DO j=jj_begin-offset,jj_end+offset990 ! DO i=ii_begin-offset,ii_end+offset991 ! ij=(j-1)*iim+i992 ! a(ij)=b(index(ij))993 !c994 RETURN995 END SUBROUTINE monGATHER996 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%997 998 subroutine monscatter(n,a,index,b)999 C1000 implicit none1001 integer N,INDEX(n),I1002 real A(n),B(n)1003 c1004 c1005 DO 100 I=1,N1006 A(INDEX(I))=B(I)1007 100 CONTINUE1008 c1009 return1010 end SUBROUTINE monscatter1011 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1012 1013 SUBROUTINE lwtr(ngrid,coef,lstrong,dup,transm)1014 IMPLICIT NONE1015 INTEGER ngrid1016 REAL coef1017 LOGICAL lstrong1018 REAL dup(ngrid),transm(ngrid)1019 INTEGER ig1020 1021 IF(lstrong) THEN1022 DO ig=1,ngrid1023 transm(ig)=exp(-coef*sqrt(dup(ig)))1024 ENDDO1025 ELSE1026 DO ig=1,ngrid1027 transm(ig)=exp(-coef*dup(ig))1028 ENDDO1029 ENDIF1030 RETURN1031 END subroutine lwtr1032 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1033 END MODULE RADIATION -
codes/icosagcm/trunk/src/surface_mod.F
r149 r186 1 MODULE SURFACE_PROCESS2 3 USE ICOSA4 USE dimphys_mod5 USE RADIATION6 DATA lmixmin,emin_turb,karman/100.,1.e-8,.4/7 8 contains9 10 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%11 12 subroutiNE vdif(ngrid,nlay,ptime,13 $ ptimestep,pcapcal,pz0,14 $ pplay,pplev,pzlay,pzlev,15 $ pu,pv,ph,ptsrf,pemis,16 $ pdufi,pdvfi,pdhfi,pfluxsrf,17 $ pdudif,pdvdif,pdhdif,pdtsrf,pq2,pq2l,18 $ lwrite)19 IMPLICIT NONE20 21 c=======================================================================22 c23 c Diffusion verticale24 c Shema implicite25 c On commence par rajouter au variables x la tendance physique26 c et on resoult en fait:27 c x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1)28 c29 c !!! attention :30 c pour utilisation sur une machine sans allocation dynamique de31 c memoires (sur SUN par exemple) il faut que ngrid soit egal32 c a ngrid.33 c34 c arguments:35 c ----------36 c37 c entree:38 c -------39 c40 c41 c=======================================================================42 43 c-----------------------------------------------------------------------44 c declarations:45 c -------------46 c47 c arguments:48 c ----------49 50 INTEGER ngrid,nlay51 REAL ptime,ptimestep52 REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)53 REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)54 REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)55 REAL ptsrf(ngrid),pemis(ngrid)56 REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)57 REAL pfluxsrf(ngrid)58 REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)59 REAL pdtsrf(ngrid),pcapcal(ngrid),pz0(ngrid)60 REAL pq2(ngrid,nlay+1),pq2l(ngrid,nlay+1)61 LOGICAL lwrite62 c63 c local:64 c ------65 66 INTEGER ilev,ig,ilay,nlev67 INTEGER unit,ierr,it1,it2,icount68 SAVE icount69 INTEGER cluvdb,putdat,putvdim,setname,setvdim70 REAL z4st,zdplanck(ngrid),zu271 REAL zkv(ngrid,nlayermx+1),zkh(ngrid,nlayermx+1)72 REAL zcdv(ngrid),zcdh(ngrid)73 REAL zu(ngrid,nlayermx),zv(ngrid,nlayermx)74 REAL zh(ngrid,nlayermx)75 REAL ztsrf2(ngrid)76 REAL z1(ngrid),z2(ngrid)77 REAL za(ngrid,nlayermx),zb(ngrid,nlayermx)78 REAL zb0(ngrid,nlayermx)79 REAL zc(ngrid,nlayermx),zd(ngrid,nlayermx)80 REAL zout_dyn(iim+1,jjm+1,nlayermx+1),zout_fi(ngrid,nlayermx+1)81 REAL zcst182 REAL karman83 84 EXTERNAL coefdifv85 EXTERNAL SSUM86 REAL SSUM87 SAVE karman88 89 DATA karman/0.4/90 DATA icount/0/91 c92 c-----------------------------------------------------------------------93 c initialisations:94 c ----------------95 96 nlev=nlay+197 98 IF(ngrid.NE.ngrid) THEN99 PRINT*,'STOP dans coefdifv'100 PRINT*,'probleme de dimensions :'101 PRINT*,'ngrid =',ngrid102 PRINT*,'ngrid =',ngrid103 STOP104 ENDIF105 106 c computation of rho*dz and dt*rho/dz=dt*rho**2 g/dp:107 c with rho=p/RT=p/ (R Theta) (p/ps)**kappa108 c ---------------------------------109 110 DO ilay=1,nlay111 DO ig=1,ngrid112 za(ig,ilay)=113 s (pplev(ig,ilay)-pplev(ig,ilay+1))/g114 ENDDO115 ENDDO116 117 zcst1=4.*g*ptimestep/(kappa*cpp)**2118 DO ilev=2,nlev-1119 DO ig=1,ngrid120 zb0(ig,ilev)=pplev(ig,ilev)*121 s (pplev(ig,1)/pplev(ig,ilev))**kappa /122 s (ph(ig,ilev-1)+ph(ig,ilev))123 zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/124 s (pplay(ig,ilev-1)-pplay(ig,ilev))125 ENDDO126 ENDDO127 DO ig=1,ngrid128 zb0(ig,1)=ptimestep*pplev(ig,1)/(kappa*cpp*ptsrf(ig))129 ENDDO130 IF(lwrite) THEN131 ig=ngrid/2+1132 PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'133 DO ilay=1,nlay134 WRITE(*,*) .01*pplay(ig,ilay),.001*pzlay(ig,ilay),135 s pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)136 ENDDO137 PRINT*,'Pression (mbar) ,altitude (km),zb'138 DO ilev=1,nlay139 WRITE(*,*) .01*pplev(ig,ilev),.001*pzlev(ig,ilev),140 s zb0(ig,ilev)141 ENDDO142 ENDIF143 144 c-----------------------------------------------------------------------145 c 2. ajout des tendances physiques:146 c ------------------------------147 148 DO ilev=1,nlay149 DO ig=1,ngrid150 zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep151 zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep152 zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep153 ENDDO154 ENDDO155 156 c-----------------------------------------------------------------------157 c 3. calcul de cd :158 c ----------------159 c160 CALL vdif_cd( ngrid,pz0,g,pzlay,pu,pv,ptsrf,ph,zcdv,zcdh)161 162 c CALL my_25(ptimestep,g,pzlev,pzlay,pu,pv,ph,zcdv,163 c a pq2,pq2l,zkv,zkh)164 165 CALL vdif_k(ngrid,nlay,166 s ptimestep,g,pzlev,pzlay,pz0,pu,pv,ph,zcdv,zkv,zkh,pq2,pq2l)167 168 DO ig=1,ngrid169 zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)170 zcdv(ig)=zcdv(ig)*sqrt(zu2)171 zcdh(ig)=zcdh(ig)*sqrt(zu2)172 ENDDO173 174 IF(lwrite) THEN175 PRINT*176 PRINT*,'Diagnostique diffusion verticale'177 PRINT*,'coefficients Cd pour v et h'178 PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)179 PRINT*,'coefficients K pour v et h'180 DO ilev=1,nlay181 PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)182 ENDDO183 ENDIF184 185 c-----------------------------------------------------------------------186 c integration verticale pour u:187 c -----------------------------188 c189 CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))190 CALL multipl(ngrid,zcdv,zb0,zb)191 DO ig=1,ngrid192 z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))193 zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig)194 zd(ig,nlay)=zb(ig,nlay)*z1(ig)195 ENDDO196 197 DO ilay=nlay-1,1,-1198 DO ig=1,ngrid199 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+200 $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))201 zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+202 $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)203 zd(ig,ilay)=zb(ig,ilay)*z1(ig)204 ENDDO205 ENDDO206 207 DO ig=1,ngrid208 zu(ig,1)=zc(ig,1)209 ENDDO210 DO ilay=2,nlay211 DO ig=1,ngrid212 zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)213 ENDDO214 ENDDO215 216 c-----------------------------------------------------------------------217 c integration verticale pour v:218 c -----------------------------219 c220 DO ig=1,ngrid221 z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))222 zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig)223 zd(ig,nlay)=zb(ig,nlay)*z1(ig)224 ENDDO225 226 DO ilay=nlay-1,1,-1227 DO ig=1,ngrid228 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+229 $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))230 zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+231 $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)232 zd(ig,ilay)=zb(ig,ilay)*z1(ig)233 ENDDO234 ENDDO235 236 DO ig=1,ngrid237 zv(ig,1)=zc(ig,1)238 ENDDO239 DO ilay=2,nlay240 DO ig=1,ngrid241 zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)242 ENDDO243 ENDDO244 245 c-----------------------------------------------------------------------246 c integration verticale pour h:247 c -----------------------------248 c249 CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))250 CALL multipl(ngrid,zcdh,zb0,zb)251 DO ig=1,ngrid252 z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))253 zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)*z1(ig)254 zd(ig,nlay)=zb(ig,nlay)*z1(ig)255 ENDDO256 257 DO ilay=nlay-1,1,-1258 DO ig=1,ngrid259 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+260 $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))261 zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+262 $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)263 zd(ig,ilay)=zb(ig,ilay)*z1(ig)264 ENDDO265 ENDDO266 267 c-----------------------------------------------------------------------268 c rajout eventuel de planck dans le shema implicite:269 c --------------------------------------------------270 271 z4st=4.*5.67e-8*ptimestep272 c z4st=0.273 DO ig=1,ngrid274 zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)275 ENDDO276 277 c-----------------------------------------------------------------------278 c calcul le l'evolution de la temperature du sol':279 c -----------------------------------------------280 281 DO ig=1,ngrid282 z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1)283 s +zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep284 z2(ig)= pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig)285 ztsrf2(ig)=z1(ig)/z2(ig)286 zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)287 pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep288 ENDDO289 290 c-----------------------------------------------------------------------291 c integration verticale finale:292 c -----------------------------293 294 DO ilay=2,nlay295 DO ig=1,ngrid296 zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zh(ig,ilay-1)297 ENDDO298 ENDDO299 300 c-----------------------------------------------------------------------301 c calcul final des tendances de la diffusion verticale:302 c -----------------------------------------------------303 304 DO ilev = 1, nlay305 DO ig=1,ngrid306 pdudif(ig,ilev)=( zu(ig,ilev)-307 $ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep308 pdvdif(ig,ilev)=( zv(ig,ilev)-309 $ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep310 pdhdif(ig,ilev)=( zh(ig,ilev)-311 $ (ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep) )/ptimestep312 ENDDO313 ENDDO314 315 IF(lwrite) THEN316 PRINT*317 PRINT*,'Diagnostique de la diffusion verticale'318 PRINT*,'h avant et apres diffusion verticale'319 PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)320 DO 3110 ilev=1,nlay321 PRINT*,ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev)322 3110 CONTINUE323 ENDIF324 c---------------------------------------------------------------------325 RETURN326 END SUBROUTINE vdif327 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%328 329 SUBROUTINE convadj(ngrid,nlay,ptimestep,330 S pplay,pplev,ppopsk,331 $ pu,pv,ph,332 $ pdufi,pdvfi,pdhfi,333 $ pduadj,pdvadj,pdhadj)334 IMPLICIT NONE335 336 c=======================================================================337 c338 c ajustement convectif sec339 c on peut ajouter les tendances pdhfi au profil pdh avant l'ajustement340 c'341 c=======================================================================342 343 c-----------------------------------------------------------------------344 c declarations:345 c -------------346 c arguments:347 c ----------348 349 INTEGER ngrid,nlay350 REAL ptimestep351 REAL ph(ngrid,nlay),pdhfi(ngrid,nlay),pdhadj(ngrid,nlay)352 REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1),ppopsk(ngrid,nlay)353 REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay)354 REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay)355 356 c local:357 c ------358 359 INTEGER ig,i,l,l1,l2,jj360 INTEGER jcnt, jadrs(ngrid)361 362 REAL*8 sig(nlayermx+1),sdsig(nlayermx),dsig(nlayermx)363 REAL*8 zu(ngrid,nlayermx),zv(ngrid,nlayermx)364 REAL*8 zh(ngrid,nlayermx)365 REAL*8 zu2(ngrid,nlayermx),zv2(ngrid,nlayermx)366 REAL*8 zh2(ngrid,nlayermx)367 REAL*8 zhm,zsm,zum,zvm,zalpha368 369 LOGICAL vtest(ngrid),down370 371 c372 c-----------------------------------------------------------------------373 c initialisation:374 c ---------------375 c376 IF(ngrid.NE.ngrid) THEN377 PRINT*378 PRINT*,'STOP dans convadj'379 PRINT*,'ngrid =',ngrid380 PRINT*,'ngrid =',ngrid381 ENDIF382 c383 c-----------------------------------------------------------------------384 c detection des profils a modifier:385 c ---------------------------------386 c si le profil est a modifier387 c (i.e. ph(niv_sup) < ph(niv_inf) )388 c alors le tableau "vtest" est mis a .TRUE. ;389 c sinon, il reste a sa valeur initiale (.FALSE.)390 c cette operation est vectorisable391 c On en profite pour copier la valeur initiale de "ph"392 c dans le champ de travail "zh"393 394 395 DO 1010 l=1,nlay396 DO 1015 ig=1,ngrid397 zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep398 zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep399 zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep400 1015 CONTINUE401 1010 CONTINUE402 403 zu2(:,:)=zu(:,:)404 zv2(:,:)=zv(:,:)405 zh2(:,:)=zh(:,:)406 407 DO 1020 ig=1,ngrid408 vtest(ig)=.FALSE.409 1020 CONTINUE410 c411 DO 1040 l=2,nlay412 DO 1060 ig=1,ngrid413 CRAY vtest(ig)=CVMGM(.TRUE. , vtest(ig),414 CRAY . zh2(ig,l)-zh2(ig,l-1))415 IF(zh2(ig,l).LT.zh2(ig,l-1)) vtest(ig)=.TRUE.416 1060 CONTINUE417 1040 CONTINUE418 c419 CRAY CALL WHENNE(ngrid, vtest, 1, 0, jadrs, jcnt)420 jcnt=0421 DO 1070 ig=1,ngrid422 IF(vtest(ig)) THEN423 jcnt=jcnt+1424 jadrs(jcnt)=ig425 ENDIF426 1070 CONTINUE427 428 429 c-----------------------------------------------------------------------430 c Ajustement des "jcnt" profils instables indices par "jadrs":431 c ------------------------------------------------------------432 c433 DO 1080 jj = 1, jcnt434 c435 i = jadrs(jj)436 c437 c Calcul des niveaux sigma sur cette colonne438 DO l=1,nlay+1439 sig(l)=pplev(i,l)/pplev(i,1)440 ENDDO441 DO l=1,nlay442 dsig(l)=sig(l)-sig(l+1)443 sdsig(l)=ppopsk(i,l)*dsig(l)444 ENDDO445 l2 = 1446 c447 c -- boucle de sondage vers le haut448 c449 cins$ Loop450 8000 CONTINUE451 c452 l2 = l2 + 1453 c454 cins$ Exit455 IF (l2 .GT. nlay) Goto 8001456 c457 IF (zh2(i, l2) .LT. zh2(i, l2-1)) THEN458 c459 c -- l2 est le niveau le plus haut de la colonne instable460 c461 l1 = l2 - 1462 l = l1463 zsm = sdsig(l2)464 zhm = zh2(i, l2)465 c466 c -- boucle de sondage vers le bas467 c468 cins$ Loop469 8020 CONTINUE470 c471 zsm = zsm + sdsig(l)472 zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm473 c474 c -- doit on etendre la colonne vers le bas ?475 c476 c_EC (M1875) 20/6/87 : AND -> AND THEN477 c478 down = .FALSE.479 IF (l1 .NE. 1) THEN !-- and then480 IF (zhm .LT. zh2(i, l1-1)) THEN481 down = .TRUE.482 END IF483 END IF484 c485 IF (down) THEN486 c487 l1 = l1 - 1488 l = l1489 c490 ELSE491 c492 c -- peut on etendre la colonne vers le haut ?493 c494 cins$ Exit495 IF (l2 .EQ. nlay) Goto 8021496 c497 cins$ Exit498 IF (zh2(i, l2+1) .GE. zhm) Goto 8021499 c500 l2 = l2 + 1501 l = l2502 c503 END IF504 c505 cins$ End Loop506 GO TO 8020507 8021 CONTINUE508 c509 c -- nouveau profil : constant (valeur moyenne)510 c511 zalpha=0.512 zum=0.513 zvm=0.514 DO 1100 l = l1, l2515 zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)516 zh2(i, l) = zhm517 zum=zum+dsig(l)*zu(i,l)518 zvm=zvm+dsig(l)*zv(i,l)519 1100 CONTINUE520 zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))521 zum=zum/(sig(l1)-sig(l2+1))522 zvm=zvm/(sig(l1)-sig(l2+1))523 IF(zalpha.GT.1.) THEN524 PRINT*,'WARNING dans convadj zalpha=',zalpha525 if(ig.eq.1) then526 print*,'Au pole nord'527 elseif (ig.eq.ngrid) then528 print*,'Au pole sud'529 else530 print*,'Point i=',531 . ig-((ig-1)/iim)*iim,'j=',(ig-1)/iim+1532 endif533 ! STOP !problem with icosa pole534 zalpha=1.535 ELSE536 c IF(zalpha.LT.0.) STOP'zalpha=0'537 IF(zalpha.LT.1.e-5) zalpha=1.e-5538 ENDIF539 DO l=l1,l2540 zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))541 zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))542 ENDDO543 544 l2 = l2 + 1545 c546 END IF547 c548 cins$ End Loop549 GO TO 8000550 8001 CONTINUE551 c552 1080 CONTINUE553 c554 DO 4000 l=1,nlay555 DO 4020 ig=1,ngrid556 pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep557 pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep558 pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep559 4020 CONTINUE560 4000 CONTINUE561 c562 RETURN563 END SUBROUTINE convadj564 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%565 566 SUBROUTINE soil(ngrid,nsoil,firstcall,ptherm_i,567 s ptimestep,ptsrf,ptsoil,568 s pcapcal,pfluxgrd)569 IMPLICIT NONE570 571 c=======================================================================572 c573 c Auteur: Frederic Hourdin 30/01/92574 c -------575 c576 c objet: computation of : the soil temperature evolution577 c ------ the surfacic heat capacity "Capcal"578 c the surface conduction flux pcapcal579 c580 c581 c Method: implicit time integration582 c -------583 c Consecutive ground temperatures are related by:584 c T(k+1) = C(k) + D(k)*T(k) (1)585 c the coefficients C and D are computed at the t-dt time-step.586 c Routine structure:587 c 1)new temperatures are computed using (1)588 c 2)C and D coefficients are computed from the new temperature589 c profile for the t+dt time-step590 c 3)the coefficients A and B are computed where the diffusive591 c fluxes at the t+dt time-step is given by592 c Fdiff = A + B Ts(t+dt)593 c or Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt594 c with F0 = A + B (Ts(t))595 c Capcal = B*dt596 c597 c Interface:598 c ----------599 c600 c Arguments:601 c ----------602 c ngird number of grid-points603 c ptimestep physical timestep (s)604 c pto(ngrid,nsoil) temperature at time-step t (K)605 c ptn(ngrid,nsoil) temperature at time step t+dt (K)606 c pcapcal(ngrid) specific heat (W*m-2*s*K-1)607 c pfluxgrd(ngrid) surface diffusive flux from ground (Wm-2)608 c609 c=======================================================================610 c declarations:611 c -------------612 c-----------------------------------------------------------------------613 c arguments614 c ---------615 616 INTEGER ngrid,nsoil617 REAL ptimestep618 REAL ptsrf(ngrid),ptsoil(ngrid,nsoilmx),ptherm_i(ngrid)619 REAL pcapcal(ngrid),pfluxgrd(ngrid)620 LOGICAL firstcall621 622 c-----------------------------------------------------------------------623 c local arrays624 c ------------625 626 INTEGER ig,jk627 REAL za(ngrid),zb(ngrid)628 REAL zdz2(nsoilmx),z1(ngrid)629 REAL min_period,dalph_soil630 631 c local saved variables:632 c ----------------------633 REAL dz1(nsoilmx),dz2(nsoilmx)634 REAL zc(ngrid,nsoilmx),zd(ngrid,nsoilmx)635 REAL lambda636 637 !!!!!!!! SARVESH !!!!!!! SAVE ATTRIBUTE638 !! SAVE dz1,dz2,zc,zd,lambda639 640 c-----------------------------------------------------------------------641 c Depthts:642 c --------643 644 REAL fz,rk,fz1,rk1,rk2645 fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)646 647 IF (firstcall) THEN648 649 c-----------------------------------------------------------------------650 c ground levels651 c grnd=z/l where l is the skin depth of the diurnal cycle:652 c --------------------------------------------------------653 654 min_period=20000.655 dalph_soil=2.656 657 OPEN(99,file='soil.def',status='old',form='formatted',err=9999)658 READ(99,*) min_period659 READ(99,*) dalph_soil660 PRINT*,'Discretization for the soil model'661 PRINT*,'First level e-folding depth',min_period,662 s ' dalph',dalph_soil663 CLOSE(99)664 9999 CONTINUE665 666 c la premiere couche represente un dixieme de cycle diurne667 fz1=sqrt(min_period/3.14)668 669 DO jk=1,nsoil670 rk1=jk671 rk2=jk-1672 dz2(jk)=fz(rk1)-fz(rk2)673 ENDDO674 DO jk=1,nsoil-1675 rk1=jk+.5676 rk2=jk-.5677 dz1(jk)=1./(fz(rk1)-fz(rk2))678 ENDDO679 lambda=fz(.5)*dz1(1)680 PRINT*,'full layers, intermediate layers (secoonds)'681 DO jk=1,nsoil682 rk=jk683 rk1=jk+.5684 rk2=jk-.5685 PRINT*,fz(rk1)*fz(rk2)*3.14,686 s fz(rk)*fz(rk)*3.14687 ENDDO688 689 c Initialisations:690 c ----------------691 692 ELSE693 c-----------------------------------------------------------------------694 c Computation of the soil temperatures using the Cgrd and Dgrd695 c coefficient computed at the previous time-step:696 c -----------------------------------------------697 698 c surface temperature699 DO ig=1,ngrid700 ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/701 s (lambda*(1.-zd(ig,1))+1.)702 ENDDO703 704 c other temperatures705 DO jk=1,nsoil-1706 DO ig=1,ngrid707 ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)708 ENDDO709 ENDDO710 711 ENDIF712 c-----------------------------------------------------------------------713 c Computation of the Cgrd and Dgrd coefficient for the next step:714 c ---------------------------------------------------------------715 716 DO jk=1,nsoil717 zdz2(jk)=dz2(jk)/ptimestep718 ENDDO719 720 DO ig=1,ngrid721 z1(ig)=zdz2(nsoil)+dz1(nsoil-1)722 zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1(ig)723 zd(ig,nsoil-1)=dz1(nsoil-1)/z1(ig)724 ENDDO725 726 DO jk=nsoil-1,2,-1727 DO ig=1,ngrid728 z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))729 zc(ig,jk-1)=730 s (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1(ig)731 zd(ig,jk-1)=dz1(jk-1)*z1(ig)732 ENDDO733 ENDDO734 735 c-----------------------------------------------------------------------736 c computation of the surface diffusive flux from ground and737 c calorific capacity of the ground:738 c ---------------------------------739 740 DO ig=1,ngrid741 pfluxgrd(ig)=ptherm_i(ig)*dz1(1)*742 s (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))743 pcapcal(ig)=ptherm_i(ig)*744 s (dz2(1)+ptimestep*(1.-zd(ig,1))*dz1(1))745 z1(ig)=lambda*(1.-zd(ig,1))+1.746 pcapcal(ig)=pcapcal(ig)/z1(ig)747 pfluxgrd(ig)=pfluxgrd(ig)748 s +pcapcal(ig)*(ptsoil(ig,1)*z1(ig)-lambda*zc(ig,1)-ptsrf(ig))749 s /ptimestep750 ENDDO751 752 RETURN753 END SUBROUTINE SOIL754 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%755 756 SUBROUTINE vdif_cd( ngrid,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh)757 IMPLICIT NONE758 c=======================================================================759 c760 c Subject: computation of the surface drag coefficient using the761 c ------- approch developed by Loui for ECMWF.762 c763 c Author: Frederic Hourdin 15 /10 /93764 c -------765 c766 c Arguments:767 c ----------768 c769 c inputs:770 c ------771 c ngrid size of the horizontal grid772 c pg gravity (m s -2)773 c pz(ngrid) height of the first atmospheric layer774 c pu(ngrid) u component of the wind in that layer775 c pv(ngrid) v component of the wind in that layer776 c pts(ngrid) surfacte temperature777 c ph(ngrid) potential temperature T*(p/ps)^kappa778 c779 c outputs:780 c --------781 c pcdv(ngrid) Cd for the wind782 c pcdh(ngrid) Cd for potential temperature783 c784 c=======================================================================785 c786 c-----------------------------------------------------------------------787 c Declarations:788 c -------------789 790 c Arguments:791 c ----------792 793 INTEGER ngrid,nlay794 REAL pz0(ngrid)795 REAL pg,pz(ngrid)796 REAL pu(ngrid),pv(ngrid)797 REAL pts(ngrid),ph(ngrid)798 REAL pcdv(ngrid),pcdh(ngrid)799 800 c Local:801 c ------802 803 INTEGER ig804 805 REAL zu2,z1,zri,zcd0,zz806 807 REAL karman,b,c,d,c2b,c3bc,c3b,z0,umin2808 LOGICAL firstcal809 DATA karman,b,c,d,umin2/.4,5.,5.,5.,1.e-12/810 DATA firstcal/.true./811 SAVE b,c,d,karman,c2b,c3bc,c3b,firstcal,umin2812 813 c-----------------------------------------------------------------------814 c couche de surface:815 c ------------------816 817 c DO ig=1,ngrid818 c zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2819 c pcdv(ig)=pz0(ig)*(1.+sqrt(zu2))820 c pcdh(ig)=pcdv(ig)821 c ENDDO822 c RETURN823 824 IF (firstcal) THEN825 c2b=2.*b826 c3bc=3.*b*c827 c3b=3.*b828 firstcal=.false.829 ENDIF830 831 c!!!! WARNING, verifier la formule originale de Louis!832 DO ig=1,ngrid833 zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2834 zri=pg*pz(ig)*(ph(ig)-pts(ig))/(ph(ig)*zu2)835 z1=1.+pz(ig)/pz0(ig)836 zcd0=karman/log(z1)837 zcd0=zcd0*zcd0*sqrt(zu2)838 IF(zri.LT.0.) THEN839 z1=b*zri/(1.+c3bc*zcd0*sqrt(-z1*zri))840 pcdv(ig)=zcd0*(1.-2.*z1)841 pcdh(ig)=zcd0*(1.-3.*z1)842 ELSE843 zz=sqrt(1.+d*zri)844 pcdv(ig)=zcd0/(1.+c2b*zri/zz)845 pcdh(ig)=zcd0/(1.+c3b*zri*zz)846 ENDIF847 ENDDO848 849 c-----------------------------------------------------------------------850 851 RETURN852 END SUBROUTINE vdif_cd853 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%854 855 SUBROUTINE vdif_k(ngrid,nlay,856 s ptimestep,pg,pzlev,pzlay,pz0,pu,pv,ph,pcdv,pkv,pkh,pq2,pq2l)857 858 IMPLICIT NONE859 860 INTEGER ngrid,nlay861 862 REAL ptimestep863 REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)864 REAL pz0(ngrid)865 REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)866 REAL pg,pcdv(ngrid)867 REAL pkv(ngrid,nlay+1),pkh(ngrid,nlay+1)868 REAL pq2(ngrid,nlay+1),pq2l(ngrid,nlay+1) !!!! SARVESH ADDED to869 870 INTEGER ig,il871 REAL zdu,zdv,zri,zdvodz2,zdz,z1,lmix872 873 REAL karman874 SAVE karman875 ! DATA lmixmin,emin_turb,karman/100.,1.e-8,.4/876 !!!!! SARVESH !!!!!!877 !Error: Host associated variable 'lmixmin' may not be in the DATA statement878 879 ! print*,'LMIXMIN',lmixmin880 DO ig=1,ngrid881 pkv(ig,1)=0.882 pkh(ig,1)=0.883 pkv(ig,nlay+1)=0.884 pkh(ig,nlay+1)=0.885 ENDDO886 c s ' zdu,zdv,zdz,zdovdz2,ph(ig,il)+ph(ig,il-1)'887 DO il=2,nlay888 DO ig=1,ngrid889 z1=pzlev(ig,il)+pz0(ig)890 lmix=karman*z1/(1.+karman*z1/lmixmin)891 c lmix=lmixmin892 c WARNING test lmix=lmixmin893 zdu=pu(ig,il)-pu(ig,il-1)894 zdv=pv(ig,il)-pv(ig,il-1)895 zdz=pzlay(ig,il)-pzlay(ig,il-1)896 zdvodz2=(zdu*zdu+zdv*zdv)/(zdz*zdz)897 IF(zdvodz2.LT.1.e-5) THEN898 pkv(ig,il)=lmix*sqrt(emin_turb)899 ELSE900 zri=2.*pg*(ph(ig,il)-ph(ig,il-1))901 s / (zdz* (ph(ig,il)+ph(ig,il-1)) *zdvodz2 )902 pkv(ig,il)=903 s lmix*sqrt(MAX(lmix*lmix*zdvodz2*(1-zri/.4),emin_turb))904 ENDIF905 pkh(ig,il)=pkv(ig,il)906 c IF(ig.EQ.ngrid/2+1) PRINT*,il,lmix,pkv(ig,il),907 c s zdu,zdv,zdz,zdvodz2,ph(ig,il)+ph(ig,il-1),908 c s lmix*lmix*zdvodz2*(1-zri/.4),emin_turb,zri,ph(ig,il)-ph(ig,il-1),909 c s ph(ig,il),ph(ig,il-1)910 ENDDO911 ENDDO912 913 RETURN914 END SUBROUTINE vdif_k915 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%916 917 SUBROUTINE multipl(n,x1,x2,y)918 IMPLICIT NONE919 c====================================================================920 c921 c multiplication de deux vecteurs922 c923 c=======================================================================924 c925 INTEGER n,i926 REAL x1(n),x2(n),y(n)927 c928 DO 10 i=1,n929 y(i)=x1(i)*x2(i)930 10 CONTINUE931 c932 RETURN933 END SUBROUTINE multipl934 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%935 END MODULE SURFACE_PROCESS -
codes/icosagcm/trunk/src/theta_rhodz.f90
r151 r186 16 16 17 17 DO ind=1,ndomain 18 IF (.NOT. assigned_domain(ind)) CYCLE 18 19 CALL swap_dimensions(ind) 19 20 CALL swap_geometry(ind) … … 39 40 40 41 DO ind=1,ndomain 42 IF (.NOT. assigned_domain(ind)) CYCLE 41 43 CALL swap_dimensions(ind) 42 44 CALL swap_geometry(ind) … … 62 64 63 65 DO ind=1,ndomain 66 IF (.NOT. assigned_domain(ind)) CYCLE 64 67 CALL swap_dimensions(ind) 65 68 CALL swap_geometry(ind) … … 82 85 INTEGER :: i,j,ij,l 83 86 REAL(rstd),ALLOCATABLE,SAVE :: p(:,:) 87 !$OMP THREADPRIVATE(p) 84 88 85 89 ALLOCATE( p(iim*jjm,llm+1)) … … 109 113 INTEGER :: i,j,ij,l 110 114 REAL(rstd),SAVE,ALLOCATABLE :: p(:,:) 115 !$OMP THREADPRIVATE(p) 111 116 112 117 ALLOCATE( p(iim*jjm,llm+1)) -
codes/icosagcm/trunk/src/time.f90
r151 r186 4 4 5 5 INTEGER,SAVE :: ncid 6 !$OMP THREADPRIVATE(ncid) 6 7 INTEGER,SAVE :: time_counter_id 8 !$OMP THREADPRIVATE(time_counter_id) 7 9 INTEGER,SAVE :: it 10 !$OMP THREADPRIVATE(it) 8 11 9 12 REAL(rstd),SAVE :: dt 13 !$OMP THREADPRIVATE(dt) 10 14 REAL(rstd),SAVE :: write_period 15 !$OMP THREADPRIVATE(write_period) 11 16 INTEGER,SAVE :: itau_out, itau_adv, itau_dissip, itau_physics, itaumax 17 !$OMP THREADPRIVATE(itau_out, itau_adv, itau_dissip, itau_physics, itaumax) 12 18 13 19 INTEGER,SAVE :: day_step,ndays 20 !$OMP THREADPRIVATE(day_step,ndays) 14 21 REAL(rstd),SAVE :: jD_ref,jH_ref 22 !$OMP THREADPRIVATE(jD_ref,jH_ref) 15 23 INTEGER,SAVE :: day_ini,day_end,annee_ref,day_ref 24 !$OMP THREADPRIVATE(day_ini,day_end,annee_ref,day_ref) 16 25 REAL(rstd),SAVE::start_time 26 !$OMP THREADPRIVATE(start_time) 17 27 CHARACTER(LEN=255) :: time_style 28 !$OMP THREADPRIVATE(time_style) 18 29 INTEGER,SAVE:: an, mois, jour 30 !$OMP THREADPRIVATE(an, mois, jour) 19 31 REAL(rstd),SAVE:: heure 32 !$OMP THREADPRIVATE(heure) 20 33 CHARACTER (LEN=10):: calend 34 !$OMP THREADPRIVATE(calend) 21 35 22 36 PUBLIC create_time_counter_header, update_time_counter, close_time_counter, init_time, & … … 32 46 SUBROUTINE init_time 33 47 USE earth_const 34 USE ioipsl48 USE getin_mod 35 49 USE mpipara 36 50 IMPLICIT NONE … … 38 52 39 53 40 54 time_style='dcmip' 41 55 CALL getin('time_style',time_style) 42 56 … … 81 95 USE netcdf_mod 82 96 USE prec 83 USE ioipsl97 USE getin_mod 84 98 USE mpipara 85 99 IMPLICIT NONE … … 88 102 REAL(rstd) :: dt 89 103 CHARACTER(LEN=255) :: time_frequency 90 104 105 CALL getin("dt",dt) 106 107 !$OMP BARRIER 108 !$OMP MASTER 91 109 IF (is_mpi_root) THEN 92 110 status = NF90_CREATE('time_counter.nc', NF90_CLOBBER, ncid) … … 104 122 status = NF90_ENDDEF(ncid) 105 123 106 CALL getin("dt",dt)107 124 status=NF90_PUT_VAR(ncid,dtid, dt) 108 125 ENDIF 109 126 it=0 127 !$OMP END MASTER 128 !$OMP BARRIER 110 129 111 130 END SUBROUTINE create_time_counter_header … … 120 139 REAL(rstd) ::time_array(1) 121 140 141 !$OMP BARRIER 122 142 !$OMP MASTER 123 143 time_array(1)=time … … 129 149 ENDIF 130 150 !$OMP END MASTER 151 !$OMP BARRIER 131 152 132 153 END SUBROUTINE update_time_counter … … 138 159 INTEGER :: status 139 160 161 !$OMP BARRIER 162 !$OMP MASTER 140 163 IF (is_mpi_root) status=NF90_CLOSE(ncid) 164 !$OMP END MASTER 165 !$OMP BARRIER 141 166 142 167 END SUBROUTINE close_time_counter -
codes/icosagcm/trunk/src/timeloop_gcm.f90
r174 r186 7 7 8 8 INTEGER, PARAMETER :: euler=1, rk4=2, mlf=3 9 INTEGER :: itau_sync=10 10 11 TYPE(t_message) :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q0 12 13 TYPE(t_field),POINTER :: f_q(:) 14 TYPE(t_field),POINTER :: f_rhodz(:), f_mass(:), f_massm1(:), f_massm2(:), f_dmass(:) 15 TYPE(t_field),POINTER :: f_phis(:), f_ps(:),f_psm1(:), f_psm2(:), f_dps(:) 16 TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:), f_du(:) 17 TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:), f_dtheta_rhodz(:) 18 TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:) 19 20 INTEGER :: nb_stage, matsuno_period, scheme 9 INTEGER, PARAMETER :: itau_sync=10 10 11 TYPE(t_message),SAVE :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q0 12 13 TYPE(t_field),POINTER,SAVE :: f_q(:) 14 TYPE(t_field),POINTER,SAVE :: f_rhodz(:), f_mass(:), f_massm1(:), f_massm2(:), f_dmass(:) 15 TYPE(t_field),POINTER,SAVE :: f_phis(:), f_ps(:),f_psm1(:), f_psm2(:), f_dps(:) 16 TYPE(t_field),POINTER,SAVE :: f_u(:),f_um1(:),f_um2(:), f_du(:) 17 TYPE(t_field),POINTER,SAVE :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:), f_dtheta_rhodz(:) 18 TYPE(t_field),POINTER,SAVE :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:) 19 20 INTEGER,SAVE :: nb_stage, matsuno_period, scheme 21 !$OMP THREADPRIVATE(nb_stage, matsuno_period, scheme) 21 22 22 23 REAL(rstd),SAVE :: jD_cur, jH_cur 24 !$OMP THREADPRIVATE(jD_cur, jH_cur) 23 25 REAL(rstd),SAVE :: start_time 24 26 !$OMP THREADPRIVATE(start_time) 25 27 CONTAINS 26 28 … … 39 41 USE transfert_mod 40 42 USE check_conserve_mod 41 USE ioipsl42 43 USE output_field_mod 43 44 USE write_field … … 47 48 48 49 !---------------------------------------------------- 49 IF (TRIM(time_style)=='lmd') Then50 51 day_step=18052 CALL getin('day_step',day_step)53 54 ndays=155 CALL getin('ndays',ndays)56 57 dt = daysec/REAL(day_step)58 itaumax = ndays*day_step59 60 calend = 'earth_360d'61 CALL getin('calend', calend)62 63 day_ini = 064 CALL getin('day_ini',day_ini)65 66 day_end = 067 CALL getin('day_end',day_end)68 69 annee_ref = 199870 CALL getin('annee_ref',annee_ref)71 72 start_time = 073 CALL getin('start_time',start_time)74 75 76 write_period=077 CALL getin('write_period',write_period)78 79 write_period=write_period/scale_factor80 itau_out=FLOOR(write_period/dt)81 82 PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out83 84 mois = 1 ; heure = 0.85 call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)86 jH_ref = jD_ref - int(jD_ref)87 jD_ref = int(jD_ref)88 89 CALL ioconf_startdate(INT(jD_ref),jH_ref)90 write(*,*)'annee_ref, mois, day_ref, heure, jD_ref'91 write(*,*)annee_ref, mois, day_ref, heure, jD_ref92 write(*,*)"ndays,day_step,itaumax,dt======>"93 write(*,*)ndays,day_step,itaumax,dt94 call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)95 write(*,*)'jD_ref+jH_ref,an, mois, jour, heure'96 write(*,*)jD_ref+jH_ref,an, mois, jour, heure97 day_end = day_ini + ndays98 END IF50 ! IF (TRIM(time_style)=='lmd') Then 51 52 ! day_step=180 53 ! CALL getin('day_step',day_step) 54 55 ! ndays=1 56 ! CALL getin('ndays',ndays) 57 58 ! dt = daysec/REAL(day_step) 59 ! itaumax = ndays*day_step 60 61 ! calend = 'earth_360d' 62 ! CALL getin('calend', calend) 63 64 ! day_ini = 0 65 ! CALL getin('day_ini',day_ini) 66 67 ! day_end = 0 68 ! CALL getin('day_end',day_end) 69 70 ! annee_ref = 1998 71 ! CALL getin('annee_ref',annee_ref) 72 73 ! start_time = 0 74 ! CALL getin('start_time',start_time) 75 76 ! 77 ! write_period=0 78 ! CALL getin('write_period',write_period) 79 ! 80 ! write_period=write_period/scale_factor 81 ! itau_out=FLOOR(write_period/dt) 82 ! 83 ! PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out 84 85 ! mois = 1 ; heure = 0. 86 ! call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref) 87 ! jH_ref = jD_ref - int(jD_ref) 88 ! jD_ref = int(jD_ref) 89 ! 90 ! CALL ioconf_startdate(INT(jD_ref),jH_ref) 91 ! write(*,*)'annee_ref, mois, day_ref, heure, jD_ref' 92 ! write(*,*)annee_ref, mois, day_ref, heure, jD_ref 93 ! write(*,*)"ndays,day_step,itaumax,dt======>" 94 ! write(*,*)ndays,day_step,itaumax,dt 95 ! call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure) 96 ! write(*,*)'jD_ref+jH_ref,an, mois, jour, heure' 97 ! write(*,*)jD_ref+jH_ref,an, mois, jour, heure 98 ! day_end = day_ini + ndays 99 ! END IF 99 100 !---------------------------------------------------- 100 101 101 102 IF (xios_output) itau_out=1 103 IF (.NOT. enable_io) itau_out=HUGE(itau_out) 102 104 103 105 ! Time-independant orography … … 178 180 CALL init_check_conserve 179 181 CALL init_physics 182 180 183 CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q) 181 184 … … 221 224 LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time 222 225 LOGICAL, PARAMETER :: check=.FALSE. 223 226 INTEGER :: start_clock 227 INTEGER :: stop_clock 228 INTEGER :: rate_clock 229 224 230 CALL caldyn_BC(f_phis, f_wflux) ! set constant values in first/last interfaces 225 231 226 ! $OMP BARRIER232 !!$OMP BARRIER 227 233 DO ind=1,ndomain 228 234 CALL swap_dimensions(ind) … … 237 243 fluxt_zero=.TRUE. 238 244 245 !$OMP MASTER 246 CALL SYSTEM_CLOCK(start_clock) 247 !$OMP END MASTER 248 249 CALL trace_on 250 239 251 DO it=0,itaumax 240 252 … … 242 254 IF (MOD(it,itau_sync)==0) THEN 243 255 CALL send_message(f_ps,req_ps0) 256 CALL wait_message(req_ps0) 244 257 CALL send_message(f_mass,req_mass0) 258 CALL wait_message(req_mass0) 245 259 CALL send_message(f_theta_rhodz,req_theta_rhodz0) 260 CALL wait_message(req_theta_rhodz0) 246 261 CALL send_message(f_u,req_u0) 262 CALL wait_message(req_u0) 247 263 CALL send_message(f_q,req_q0) 248 CALL wait_message(req_ps0)249 CALL wait_message(req_mass0)250 CALL wait_message(req_theta_rhodz0)251 CALL wait_message(req_u0)252 264 CALL wait_message(req_q0) 265 266 ! CALL wait_message(req_ps0) 267 ! CALL wait_message(req_mass0) 268 ! CALL wait_message(req_theta_rhodz0) 269 ! CALL wait_message(req_u0) 270 ! CALL wait_message(req_q0) 253 271 ENDIF 254 255 ! IF (is_mpi_root) PRINT *,"It No :",It," t :",dt*It 272 273 !$OMP MASTER 274 IF (is_mpi_root) PRINT *,"It No :",It," t :",dt*It 275 !$OMP END MASTER 256 276 IF (mod(it,itau_out)==0 ) THEN 257 277 CALL update_time_counter(dt*it) … … 279 299 280 300 IF (MOD(it+1,itau_dissip)==0) THEN 301 ! CALL send_message(f_ps,req_ps) 302 ! CALL wait_message(req_ps) 303 281 304 IF(caldyn_eta==eta_mass) THEN 282 305 DO ind=1,ndomain 306 IF (.NOT. assigned_domain(ind)) CYCLE 283 307 CALL swap_dimensions(ind) 284 308 CALL swap_geometry(ind) … … 287 311 END DO 288 312 ENDIF 313 ! CALL send_message(f_mass,req_mass) 314 ! CALL wait_message(req_mass) 289 315 CALL dissip(f_u,f_du,f_mass,f_phis, f_theta_rhodz,f_dtheta_rhodz) 316 ! CALL send_message(f_mass,req_mass) 317 ! CALL wait_message(req_mass) 290 318 CALL euler_scheme(.FALSE.) ! update only u, theta 291 319 END IF … … 299 327 IF (check) THEN 300 328 DO ind=1,ndomain 329 IF (.NOT. assigned_domain(ind)) CYCLE 301 330 CALL swap_dimensions(ind) 302 331 CALL swap_geometry(ind) … … 316 345 ! jH_cur = jH_cur - int(jH_cur) 317 346 CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q) 347 318 348 ENDDO 319 349 350 CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 351 352 !$OMP MASTER 353 CALL SYSTEM_CLOCK(stop_clock) 354 CALL SYSTEM_CLOCK(count_rate=rate_clock) 355 356 IF (mpi_rank==0) THEN 357 PRINT *,"Time elapsed : ",(stop_clock-start_clock)*1./rate_clock 358 ENDIF 359 !$OMP END MASTER 320 360 321 361 CONTAINS … … 329 369 330 370 DO ind=1,ndomain 371 IF (.NOT. assigned_domain(ind)) CYCLE 331 372 CALL swap_dimensions(ind) 332 373 CALL swap_geometry(ind) … … 390 431 391 432 DO ind=1,ndomain 433 IF (.NOT. assigned_domain(ind)) CYCLE 392 434 CALL swap_dimensions(ind) 393 435 CALL swap_geometry(ind) … … 408 450 ENDDO 409 451 ENDIF 410 CALL send_message(f_ps,req_ps) 452 ! CALL send_message(f_ps,req_ps) 453 !ym no overlap for now 454 ! CALL wait_message(req_ps) 411 455 412 456 ELSE ! Lagrangian coordinate, deal with mass 413 457 DO ind=1,ndomain 458 IF (.NOT. assigned_domain(ind)) CYCLE 414 459 CALL swap_dimensions(ind) 415 460 CALL swap_geometry(ind) … … 433 478 ENDDO 434 479 END DO 435 CALL send_message(f_mass,req_mass) 480 ! CALL send_message(f_mass,req_mass) 481 !ym no overlap for now 482 ! CALL wait_message(req_mass) 436 483 437 484 END IF … … 439 486 ! now deal with other prognostic variables 440 487 DO ind=1,ndomain 488 IF (.NOT. assigned_domain(ind)) CYCLE 441 489 CALL swap_dimensions(ind) 442 490 CALL swap_geometry(ind) … … 488 536 tau = dt/nb_stage 489 537 DO ind=1,ndomain 538 IF (.NOT. assigned_domain(ind)) CYCLE 490 539 CALL swap_dimensions(ind) 491 540 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/trace.F90
r151 r186 17 17 END SUBROUTINE init_trace 18 18 19 SUBROUTINE trace_on 20 IMPLICIT NONE 21 #ifdef VTRACE 22 #include <vt_user.inc> 23 #endif 24 !$OMP MASTER 25 #ifdef VTRACE 26 VT_ON() 27 #endif 28 !$OMP END MASTER 29 30 END SUBROUTINE trace_on 19 31 32 SUBROUTINE trace_off 33 IMPLICIT NONE 34 #ifdef VTRACE 35 #include <vt_user.inc> 36 #endif 37 38 !$OMP MASTER 39 #ifdef VTRACE 40 VT_OFF() 41 #endif 42 !$OMP END MASTER 43 44 END SUBROUTINE trace_off 45 20 46 SUBROUTINE trace_start(name) 21 47 IMPLICIT NONE … … 49 75 END SUBROUTINE trace_end 50 76 77 SUBROUTINE trace_start2(name) 78 IMPLICIT NONE 79 CHARACTER(LEN=*),INTENT(IN) :: name 80 #ifdef VTRACE 81 #include <vt_user.inc> 82 #endif 83 84 !$OMP MASTER 85 #ifdef VTRACE 86 VT_USER_START(name) 87 #endif 88 !$OMP END MASTER 89 90 END SUBROUTINE trace_start2 91 92 SUBROUTINE trace_end2(name) 93 IMPLICIT NONE 94 #ifdef VTRACE 95 #include <vt_user.inc> 96 #endif 97 98 CHARACTER(LEN=*),INTENT(IN) :: name 99 100 !$OMP MASTER 101 #ifdef VTRACE 102 VT_USER_END(name) 103 #endif 104 !$OMP END MASTER 105 106 END SUBROUTINE trace_end2 107 108 109 51 110 SUBROUTINE Marker(name) 52 111 IMPLICIT NONE -
codes/icosagcm/trunk/src/transfert_mpi.f90
r176 r186 8 8 INTEGER :: domain 9 9 INTEGER :: rank 10 INTEGER :: tag 10 11 INTEGER :: size 12 INTEGER :: offset 13 INTEGER :: ireq 11 14 INTEGER,POINTER :: buffer(:) 12 REAL,POINTER :: buffer_r2(:) 13 REAL,POINTER :: buffer_r3(:,:) 14 REAL,POINTER :: buffer_r4(:,:,:) 15 REAL,POINTER :: buffer_r(:) 15 16 INTEGER,POINTER :: src_value(:) 16 17 END TYPE array 17 18 18 19 TYPE t_buffer 19 REAL,POINTER :: r 2(:)20 REAL,POINTER :: r3(:,:)21 REAL,POINTER :: r4(:,:,:)20 REAL,POINTER :: r(:) 21 INTEGER :: size 22 INTEGER :: rank 22 23 END TYPE t_buffer 23 24 … … 39 40 INTEGER :: nsend 40 41 INTEGER :: nreq_mpi 42 INTEGER :: nreq_send 43 INTEGER :: nreq_recv 41 44 TYPE(ARRAY),POINTER :: send(:) 42 45 END TYPE t_request 43 46 44 TYPE(t_request),POINTER :: req_i1(:) 45 TYPE(t_request),POINTER :: req_e1_scal(:) 46 TYPE(t_request),POINTER :: req_e1_vect(:) 47 48 TYPE(t_request),POINTER :: req_i0(:) 49 TYPE(t_request),POINTER :: req_e0_scal(:) 50 TYPE(t_request),POINTER :: req_e0_vect(:) 47 TYPE(t_request),SAVE,POINTER :: req_i1(:) 48 TYPE(t_request),SAVE,POINTER :: req_e1_scal(:) 49 TYPE(t_request),SAVE,POINTER :: req_e1_vect(:) 50 51 TYPE(t_request),SAVE,POINTER :: req_i0(:) 52 TYPE(t_request),SAVE,POINTER :: req_e0_scal(:) 53 TYPE(t_request),SAVE,POINTER :: req_e0_vect(:) 54 55 TYPE t_reorder 56 INTEGER :: ind 57 INTEGER :: rank 58 INTEGER :: tag 59 INTEGER :: isend 60 END TYPE t_reorder 51 61 52 62 TYPE t_message 53 63 TYPE(t_request), POINTER :: request(:) 54 64 INTEGER :: nreq 65 INTEGER :: nreq_send 66 INTEGER :: nreq_recv 67 TYPE(t_reorder), POINTER :: reorder(:) 55 68 INTEGER, POINTER :: mpi_req(:) 56 69 INTEGER, POINTER :: status(:,:) … … 62 75 END TYPE t_message 63 76 64 INTEGER,SAVE :: message_number=0 ;65 66 77 CONTAINS 67 78 79 68 80 SUBROUTINE init_transfert 69 81 USE domain_mod … … 72 84 USE metric 73 85 USE mpipara 86 USE mpi_mod 74 87 IMPLICIT NONE 75 88 INTEGER :: ind,i,j 89 LOGICAL ::ok 76 90 77 91 CALL create_request(field_t,req_i1) … … 410 424 INTEGER :: nb_domain_recv(0:mpi_size-1) 411 425 INTEGER :: nb_domain_send(0:mpi_size-1) 426 INTEGER :: tag_rank(0:mpi_size-1) 412 427 INTEGER :: nb_data_domain_recv(ndomain_glo) 413 428 INTEGER :: list_domain_recv(ndomain_glo) … … 415 430 INTEGER :: list_domain(ndomain) 416 431 417 INTEGER :: rank,i,j 432 INTEGER :: rank,i,j,pos 418 433 INTEGER :: size_,ind_glo,ind_loc, ind_src 419 INTEGER :: isend, irecv, ireq, nreq 434 INTEGER :: isend, irecv, ireq, nreq, nsend, nrecv 420 435 INTEGER, ALLOCATABLE :: mpi_req(:) 421 436 INTEGER, ALLOCATABLE :: status(:,:) 422 437 INTEGER, ALLOCATABLE :: rank_list(:) 438 INTEGER, ALLOCATABLE :: offset(:) 439 LOGICAL,PARAMETER :: debug = .FALSE. 440 441 423 442 IF (.NOT. using_mpi) RETURN 424 443 … … 428 447 nb_data_domain_recv(:) = 0 429 448 nb_domain_recv(:) = 0 449 tag_rank(:)=0 430 450 431 451 DO i=1,req%size … … 486 506 ALLOCATE(status(MPI_STATUS_SIZE,nreq)) 487 507 508 488 509 ireq=0 489 510 DO ind_loc=1,ndomain … … 492 513 ireq=ireq+1 493 514 CALL MPI_ISEND(req%recv(irecv)%domain,1,MPI_INTEGER,req%recv(irecv)%rank,0,comm_icosa, mpi_req(ireq),ierr) 515 IF (debug) PRINT *,"Isend ",req%recv(irecv)%domain, "from ", mpi_rank, "to ",req%recv(irecv)%rank, "tag ",0 494 516 ENDDO 495 517 ENDDO 496 518 519 IF (debug) PRINT *,"------------" 497 520 j=0 498 521 DO rank=0,mpi_size-1 … … 501 524 ireq=ireq+1 502 525 CALL MPI_IRECV(list_domain_send(j),1,MPI_INTEGER,rank,0,comm_icosa, mpi_req(ireq),ierr) 526 IF (debug) PRINT *,"IRecv ",list_domain_send(j), "from ", rank, "to ",mpi_rank, "tag ",0 503 527 ENDDO 504 528 ENDDO 529 IF (debug) PRINT *,"------------" 505 530 506 531 CALL MPI_WAITALL(nreq,mpi_req,status,ierr) … … 517 542 ALLOCATE(req%send(req%nsend)) 518 543 ENDDO 544 545 IF (debug) PRINT *,"------------" 519 546 520 547 ireq=0 … … 525 552 ireq=ireq+1 526 553 CALL MPI_ISEND(mpi_rank,1,MPI_INTEGER,req%recv(irecv)%rank,req%recv(irecv)%domain,comm_icosa, mpi_req(ireq),ierr) 554 IF (debug) PRINT *,"Isend ",mpi_rank, "from ", mpi_rank, "to ",req%recv(irecv)%rank,"tag ",req%recv(irecv)%domain 527 555 ENDDO 556 IF (debug) PRINT *,"------------" 528 557 529 558 DO isend=1,req%nsend 530 559 ireq=ireq+1 531 560 CALL MPI_IRECV(req%send(isend)%rank,1,MPI_INTEGER,MPI_ANY_SOURCE,ind_loc,comm_icosa, mpi_req(ireq),ierr) 561 IF (debug) PRINT *,"IRecv ",req%send(isend)%rank, "from ", "xxx", "to ",mpi_rank, "tag ",ind_loc 532 562 ENDDO 533 563 ENDDO 534 564 565 IF (debug) PRINT *,"------------" 566 535 567 CALL MPI_WAITALL(nreq,mpi_req,status,ierr) 536 568 CALL MPI_BARRIER(comm_icosa,ierr) 569 570 IF (debug) PRINT *,"------------" 537 571 538 572 ireq=0 … … 543 577 ireq=ireq+1 544 578 CALL MPI_ISEND(ind_loc,1,MPI_INTEGER,req%recv(irecv)%rank,req%recv(irecv)%domain,comm_icosa, mpi_req(ireq),ierr) 579 IF (debug) PRINT *,"Isend ",ind_loc, "from ", mpi_rank, "to ",req%recv(irecv)%rank,"tag ",req%recv(irecv)%domain 545 580 ENDDO 581 582 IF (debug) PRINT *,"------------" 546 583 547 584 DO isend=1,req%nsend 548 585 ireq=ireq+1 549 586 CALL MPI_IRECV(req%send(isend)%domain,1,MPI_INTEGER,req%send(isend)%rank,ind_loc,comm_icosa, mpi_req(ireq),ierr) 587 IF (debug) PRINT *,"IRecv ",req%send(isend)%domain, "from ", req%send(isend)%rank, "to ",mpi_rank, "tag ",ind_loc 550 588 ENDDO 551 589 ENDDO 590 IF (debug) PRINT *,"------------" 552 591 553 592 CALL MPI_WAITALL(nreq,mpi_req,status,ierr) 554 593 CALL MPI_BARRIER(comm_icosa,ierr) 594 IF (debug) PRINT *,"------------" 595 596 ireq=0 597 DO ind_loc=1,ndomain 598 req=>request(ind_loc) 599 600 DO irecv=1,req%nrecv 601 ireq=ireq+1 602 req%recv(irecv)%tag=tag_rank(req%recv(irecv)%rank) 603 tag_rank(req%recv(irecv)%rank)=tag_rank(req%recv(irecv)%rank)+1 604 CALL MPI_ISEND(req%recv(irecv)%tag,1,MPI_INTEGER,req%recv(irecv)%rank,req%recv(irecv)%domain,comm_icosa, mpi_req(ireq),ierr) 605 IF (debug) PRINT *,"Isend ",req%recv(irecv)%tag, "from ", mpi_rank, "to ",req%recv(irecv)%rank,"tag ",req%recv(irecv)%domain 606 ENDDO 607 IF (debug) PRINT *,"------------" 608 609 DO isend=1,req%nsend 610 ireq=ireq+1 611 CALL MPI_IRECV(req%send(isend)%tag,1,MPI_INTEGER,req%send(isend)%rank,ind_loc,comm_icosa, mpi_req(ireq),ierr) 612 IF (debug) PRINT *,"IRecv ",req%send(isend)%tag, "from ", req%send(isend)%rank, "to ",mpi_rank, "tag ",ind_loc 613 ENDDO 614 ENDDO 615 IF (debug) PRINT *,"------------" 616 617 CALL MPI_WAITALL(nreq,mpi_req,status,ierr) 618 CALL MPI_BARRIER(comm_icosa,ierr) 619 620 621 IF (debug) PRINT *,"------------" 555 622 556 623 ireq=0 … … 560 627 DO irecv=1,req%nrecv 561 628 ireq=ireq+1 562 CALL MPI_ISEND(req%recv(irecv)%size,1,MPI_INTEGER,req%recv(irecv)%rank,req%recv(irecv)%domain,comm_icosa, mpi_req(ireq),ierr) 629 CALL MPI_ISEND(req%recv(irecv)%size,1,MPI_INTEGER,req%recv(irecv)%rank,req%recv(irecv)%tag,comm_icosa, mpi_req(ireq),ierr) 630 IF (debug) PRINT *,"Isend ",req%recv(irecv)%size, "from ", mpi_rank, "to ",req%recv(irecv)%rank,"tag ",req%recv(irecv)%domain 563 631 ENDDO 632 IF (debug) PRINT *,"------------" 564 633 565 634 DO isend=1,req%nsend 566 635 ireq=ireq+1 567 CALL MPI_IRECV(req%send(isend)%size,1,MPI_INTEGER,req%send(isend)%rank,ind_loc,comm_icosa, mpi_req(ireq),ierr) 636 CALL MPI_IRECV(req%send(isend)%size,1,MPI_INTEGER,req%send(isend)%rank,req%send(isend)%tag,comm_icosa, mpi_req(ireq),ierr) 637 IF (debug) PRINT *,"IRecv ",req%send(isend)%size, "from ", req%send(isend)%rank, "to ",mpi_rank, "tag ",ind_loc 568 638 ENDDO 569 639 ENDDO … … 578 648 ireq=ireq+1 579 649 CALL MPI_ISEND(req%recv(irecv)%value,req%recv(irecv)%size,MPI_INTEGER,& 580 req%recv(irecv)%rank,req%recv(irecv)% domain,comm_icosa, mpi_req(ireq),ierr)650 req%recv(irecv)%rank,req%recv(irecv)%tag,comm_icosa, mpi_req(ireq),ierr) 581 651 ENDDO 582 652 … … 585 655 ALLOCATE(req%send(isend)%value(req%send(isend)%size)) 586 656 CALL MPI_IRECV(req%send(isend)%value,req%send(isend)%size,MPI_INTEGER,& 587 req%send(isend)%rank, ind_loc,comm_icosa, mpi_req(ireq),ierr)657 req%send(isend)%rank,req%send(isend)%tag,comm_icosa, mpi_req(ireq),ierr) 588 658 ENDDO 589 659 ENDDO … … 600 670 ENDDO 601 671 ENDDO 602 603 ! domain is on the same mpi process 672 673 674 ! domain is on the same mpi process => copie memory to memory 604 675 605 676 DO ind_loc=1,ndomain … … 611 682 req_src=>request(req%recv(irecv)%domain) 612 683 DO isend=1,req_src%nsend 613 IF (req_src%send(isend)%rank==mpi_rank .AND. req_src%send(isend)% domain==ind_loc) THEN684 IF (req_src%send(isend)%rank==mpi_rank .AND. req_src%send(isend)%tag==req%recv(irecv)%tag) THEN 614 685 req%recv(irecv)%src_value => req_src%send(isend)%value 615 686 IF ( size(req%recv(irecv)%value) /= size(req_src%send(isend)%value)) THEN 687 PRINT *,ind_loc, irecv, isend, size(req%recv(irecv)%value), size(req_src%send(isend)%value) 616 688 STOP "size(req%recv(irecv)%value) /= size(req_src%send(isend)%value" 617 689 ENDIF … … 624 696 625 697 ! true number of mpi request 698 699 request(:)%nreq_mpi=0 700 request(:)%nreq_send=0 701 request(:)%nreq_recv=0 702 ALLOCATE(rank_list(sum(request(:)%nsend))) 703 ALLOCATE(offset(sum(request(:)%nsend))) 704 offset(:)=0 705 706 nsend=0 626 707 DO ind_loc=1,ndomain 627 708 req=>request(ind_loc) 628 req%nreq_mpi=0629 709 630 710 DO isend=1,req%nsend 631 IF (req%send(isend)%rank/=mpi_rank .OR. .TRUE.) req%nreq_mpi=req%nreq_mpi+1 711 IF (req%send(isend)%rank/=mpi_rank) THEN 712 pos=0 713 DO i=1,nsend 714 IF (req%send(isend)%rank==rank_list(i)) EXIT 715 pos=pos+1 716 ENDDO 717 718 IF (pos==nsend) THEN 719 nsend=nsend+1 720 req%nreq_mpi=req%nreq_mpi+1 721 req%nreq_send=req%nreq_send+1 722 IF (mpi_threading_mode==MPI_THREAD_FUNNELED) THEN 723 rank_list(nsend)=req%send(isend)%rank 724 ELSE 725 rank_list(nsend)=-1 726 ENDIF 727 ENDIF 728 729 pos=pos+1 730 req%send(isend)%ireq=pos 731 req%send(isend)%offset=offset(pos) 732 offset(pos)=offset(pos)+req%send(isend)%size 733 ENDIF 632 734 ENDDO 633 735 ENDDO 736 737 DEALLOCATE(rank_list) 738 DEALLOCATE(offset) 739 740 ALLOCATE(rank_list(sum(request(:)%nrecv))) 741 ALLOCATE(offset(sum(request(:)%nrecv))) 742 offset(:)=0 743 744 nrecv=0 745 DO ind_loc=1,ndomain 746 req=>request(ind_loc) 747 634 748 DO irecv=1,req%nrecv 635 IF (req%recv(irecv)%rank/=mpi_rank .OR. .TRUE.) req%nreq_mpi=req%nreq_mpi+1 749 IF (req%recv(irecv)%rank/=mpi_rank) THEN 750 pos=0 751 DO i=1,nrecv 752 IF (req%recv(irecv)%rank==rank_list(i)) EXIT 753 pos=pos+1 754 ENDDO 755 756 IF (pos==nrecv) THEN 757 nrecv=nrecv+1 758 req%nreq_mpi=req%nreq_mpi+1 759 req%nreq_recv=req%nreq_recv+1 760 IF (mpi_threading_mode==MPI_THREAD_FUNNELED) THEN 761 rank_list(nrecv)=req%recv(irecv)%rank 762 ELSE 763 rank_list(nrecv)=-1 764 ENDIF 765 ENDIF 766 767 pos=pos+1 768 req%recv(irecv)%ireq=nsend+pos 769 req%recv(irecv)%offset=offset(pos) 770 offset(pos)=offset(pos)+req%recv(irecv)%size 771 ENDIF 636 772 ENDDO 637 638 773 ENDDO 774 775 ! get the offsets 776 777 ireq=0 778 DO ind_loc=1,ndomain 779 req=>request(ind_loc) 780 781 DO irecv=1,req%nrecv 782 ireq=ireq+1 783 CALL MPI_ISEND(req%recv(irecv)%offset,1,MPI_INTEGER,& 784 req%recv(irecv)%rank,req%recv(irecv)%tag,comm_icosa, mpi_req(ireq),ierr) 785 ENDDO 786 787 DO isend=1,req%nsend 788 ireq=ireq+1 789 CALL MPI_IRECV(req%send(isend)%offset,1,MPI_INTEGER,& 790 req%send(isend)%rank,req%send(isend)%tag,comm_icosa, mpi_req(ireq),ierr) 791 ENDDO 792 ENDDO 793 794 CALL MPI_WAITALL(nreq,mpi_req,status,ierr) 795 639 796 640 797 END SUBROUTINE Finalize_request … … 701 858 END SUBROUTINE transfert_message_seq 702 859 860 861 703 862 704 863 SUBROUTINE init_message_mpi(field,request, message) … … 717 876 TYPE(t_request),POINTER :: req 718 877 INTEGER :: irecv,isend 719 INTEGER :: ireq,nreq 878 INTEGER :: ireq,nreq, nreq_send 720 879 INTEGER :: ind 721 880 INTEGER :: dim3,dim4 722 881 INTEGER :: i,j 882 INTEGER :: message_number 883 ! TYPE(t_reorder),POINTER :: reorder(:) 884 ! TYPE(t_reorder) :: reorder_swap 885 886 !$OMP BARRIER 723 887 !$OMP MASTER 724 888 message%number=message_number 725 889 message_number=message_number+1 726 890 IF (message_number==100) message_number=0 727 891 892 728 893 message%request=>request 729 nreq=sum(request(:)%nsend)+sum(request(:)%nrecv)730 ! message%nreq=nreq731 894 message%nreq=sum(message%request(:)%nreq_mpi) 895 message%nreq_send=sum(message%request(:)%nreq_send) 896 message%nreq_recv=sum(message%request(:)%nreq_recv) 897 nreq=message%nreq 898 732 899 ALLOCATE(message%mpi_req(nreq)) 733 900 ALLOCATE(message%buffers(nreq)) 734 901 ALLOCATE(message%status(MPI_STATUS_SIZE,nreq)) 735 902 message%buffers(:)%size=0 736 903 message%pending=.FALSE. 737 904 message%completed=.FALSE. 738 905 906 DO ind=1,ndomain 907 req=>request(ind) 908 DO isend=1,req%nsend 909 IF (req%send(isend)%rank/=mpi_rank) THEN 910 ireq=req%send(isend)%ireq 911 message%buffers(ireq)%size=message%buffers(ireq)%size+req%send(isend)%size 912 message%buffers(ireq)%rank=req%send(isend)%rank 913 ENDIF 914 ENDDO 915 DO irecv=1,req%nrecv 916 IF (req%recv(irecv)%rank/=mpi_rank) THEN 917 ireq=req%recv(irecv)%ireq 918 message%buffers(ireq)%size=message%buffers(ireq)%size+req%recv(irecv)%size 919 message%buffers(ireq)%rank=req%recv(irecv)%rank 920 ENDIF 921 ENDDO 922 ENDDO 923 924 739 925 IF (field(1)%data_type==type_real) THEN 740 926 741 927 IF (field(1)%ndim==2) THEN 742 743 ireq=0 744 DO ind=1,ndomain 745 req=>request(ind) 746 747 DO isend=1,req%nsend 748 ireq=ireq+1 749 send=>req%send(isend) 750 CALL allocate_mpi_buffer(message%buffers(ireq)%r2,send%size) 751 ENDDO 752 753 DO irecv=1,req%nrecv 754 ireq=ireq+1 755 recv=>req%recv(irecv) 756 CALL allocate_mpi_buffer(message%buffers(ireq)%r2,recv%size) 757 ENDDO 758 759 ENDDO 760 928 929 DO ireq=1,message%nreq 930 CALL allocate_mpi_buffer(message%buffers(ireq)%r,message%buffers(ireq)%size) 931 ENDDO 761 932 762 933 ELSE IF (field(1)%ndim==3) THEN 763 764 ireq=0 765 DO ind=1,ndomain 766 dim3=size(field(ind)%rval3d,2) 767 req=>request(ind) 768 769 DO isend=1,req%nsend 770 ireq=ireq+1 771 send=>req%send(isend) 772 CALL allocate_mpi_buffer(message%buffers(ireq)%r3,send%size,dim3) 773 ENDDO 774 775 DO irecv=1,req%nrecv 776 ireq=ireq+1 777 recv=>req%recv(irecv) 778 CALL allocate_mpi_buffer(message%buffers(ireq)%r3,recv%size,dim3) 779 780 ENDDO 781 782 ENDDO 783 784 934 935 dim3=size(field(1)%rval3d,2) 936 DO ireq=1,message%nreq 937 message%buffers(ireq)%size=message%buffers(ireq)%size*dim3 938 CALL allocate_mpi_buffer(message%buffers(ireq)%r,message%buffers(ireq)%size) 939 ENDDO 940 785 941 ELSE IF (field(1)%ndim==4) THEN 786 787 ireq=0 788 DO ind=1,ndomain 789 dim3=size(field(ind)%rval4d,2) 790 dim4=size(field(ind)%rval4d,3) 791 req=>request(ind) 792 793 DO isend=1,req%nsend 794 ireq=ireq+1 795 send=>req%send(isend) 796 CALL allocate_mpi_buffer(message%buffers(ireq)%r4,send%size,dim3,dim4) 797 ENDDO 798 799 DO irecv=1,req%nrecv 800 ireq=ireq+1 801 recv=>req%recv(irecv) 802 CALL allocate_mpi_buffer(message%buffers(ireq)%r4,recv%size,dim3,dim4) 803 ENDDO 804 805 ENDDO 806 942 dim3=size(field(1)%rval4d,2) 943 dim4=size(field(1)%rval4d,3) 944 DO ireq=1,message%nreq 945 message%buffers(ireq)%size=message%buffers(ireq)%size*dim3*dim4 946 CALL allocate_mpi_buffer(message%buffers(ireq)%r,message%buffers(ireq)%size) 947 ENDDO 807 948 ENDIF 808 949 ENDIF 950 951 952 953 ! ! Reorder the request, so recv request are done in the same order than send request 954 955 ! nreq_send=sum(request(:)%nsend) 956 ! message%nreq_send=nreq_send 957 ! ALLOCATE(message%reorder(nreq_send)) 958 ! reorder=>message%reorder 959 ! ireq=0 960 ! DO ind=1,ndomain 961 ! req=>request(ind) 962 ! DO isend=1,req%nsend 963 ! ireq=ireq+1 964 ! reorder(ireq)%ind=ind 965 ! reorder(ireq)%isend=isend 966 ! reorder(ireq)%tag=req%send(isend)%tag 967 ! ENDDO 968 ! ENDDO 969 970 ! ! do a very very bad sort 971 ! DO i=1,nreq_send-1 972 ! DO j=i+1,nreq_send 973 ! IF (reorder(i)%tag > reorder(j)%tag) THEN 974 ! reorder_swap=reorder(i) 975 ! reorder(i)=reorder(j) 976 ! reorder(j)=reorder_swap 977 ! ENDIF 978 ! ENDDO 979 ! ENDDO 980 ! PRINT *,"reorder ",reorder(:)%tag 981 982 809 983 !$OMP END MASTER 810 984 !$OMP BARRIER 985 811 986 END SUBROUTINE init_message_mpi 987 988 SUBROUTINE Finalize_message_mpi(field,message) 989 USE field_mod 990 USE domain_mod 991 USE mpi_mod 992 USE mpipara 993 USE mpi_mod 994 IMPLICIT NONE 995 TYPE(t_field),POINTER :: field(:) 996 TYPE(t_message) :: message 997 998 TYPE(t_request),POINTER :: req 999 INTEGER :: irecv,isend 1000 INTEGER :: ireq,nreq 1001 INTEGER :: ind 1002 1003 !$OMP BARRIER 1004 !$OMP MASTER 1005 1006 1007 IF (field(1)%data_type==type_real) THEN 1008 1009 IF (field(1)%ndim==2) THEN 1010 1011 DO ireq=1,message%nreq 1012 CALL free_mpi_buffer(message%buffers(ireq)%r) 1013 ENDDO 1014 1015 ELSE IF (field(1)%ndim==3) THEN 1016 1017 DO ireq=1,message%nreq 1018 CALL free_mpi_buffer(message%buffers(ireq)%r) 1019 ENDDO 1020 1021 ELSE IF (field(1)%ndim==4) THEN 1022 1023 DO ireq=1,message%nreq 1024 CALL free_mpi_buffer(message%buffers(ireq)%r) 1025 ENDDO 1026 1027 ENDIF 1028 ENDIF 1029 1030 1031 1032 !$OMP END MASTER 1033 !$OMP BARRIER 1034 1035 1036 END SUBROUTINE Finalize_message_mpi 1037 1038 812 1039 813 1040 SUBROUTINE barrier … … 845 1072 REAL(rstd),POINTER :: rval3d(:,:), src_rval3d(:,:) 846 1073 REAL(rstd),POINTER :: rval4d(:,:,:), src_rval4d(:,:,:) 847 REAL(rstd),POINTER :: buffer_r2(:) 848 REAL(rstd),POINTER :: buffer_r3(:,:) 849 REAL(rstd),POINTER :: buffer_r4(:,:,:) 850 INTEGER,POINTER :: value(:) 851 INTEGER,POINTER :: sgn(:) 852 TYPE(ARRAY),POINTER :: recv,send 853 TYPE(t_request),POINTER :: req 854 INTEGER, ALLOCATABLE :: mpi_req(:) 855 INTEGER, ALLOCATABLE :: status(:,:) 856 INTEGER :: irecv,isend 857 INTEGER :: ireq,ireq_mpi,nreq 858 INTEGER :: ind,n,l,m 859 INTEGER :: dim3,dim4 860 INTEGER,POINTER :: src_value(:) 861 INTEGER,POINTER :: sign(:) 862 863 !$OMP BARRIER 864 865 CALL trace_start("transfert_mpi") 866 867 ! nreq=message%nreq 868 message%field=>field 869 870 !$OMP MASTER 871 IF (message%nreq>0) THEN 872 message%completed=.FALSE. 873 message%pending=.TRUE. 874 ELSE 875 message%completed=.TRUE. 876 message%pending=.FALSE. 877 ENDIF 878 879 !$OMP END MASTER 880 881 IF (field(1)%data_type==type_real) THEN 882 IF (field(1)%ndim==2) THEN 883 884 ireq=0 885 ireq_mpi=0 886 DO ind=1,ndomain 887 rval2d=>field(ind)%rval2d 888 889 req=>message%request(ind) 890 DO isend=1,req%nsend 891 ireq=ireq+1 892 send=>req%send(isend) 893 value=>send%value 894 895 896 IF (send%rank/=mpi_rank .OR. .TRUE.) THEN 897 ireq_mpi=ireq_mpi+1 898 buffer_r2=>message%buffers(ireq)%r2 899 CALL trace_in 900 901 !$OMP DO SCHEDULE(STATIC) 902 DO n=1,send%size 903 buffer_r2(n)=rval2d(value(n)) 904 ENDDO 905 906 CALL trace_out 907 908 !$OMP MASTER 909 CALL MPI_ISSEND(buffer_r2,send%size,MPI_REAL8,send%rank,ind+100*message%number,comm_icosa, message%mpi_req(ireq_mpi),ierr) 910 !$OMP END MASTER 911 912 ENDIF 913 ENDDO 914 915 DO irecv=1,req%nrecv 916 ireq=ireq+1 917 recv=>req%recv(irecv) 918 919 IF (recv%rank==mpi_rank .AND. .FALSE.) THEN 920 value=>recv%value 921 src_value => recv%src_value 922 src_rval2d=>field(recv%domain)%rval2d 923 sgn=>recv%sign 924 !$OMP DO SCHEDULE(STATIC) 925 DO n=1,recv%size 926 rval2d(value(n))=src_rval2d(src_value(n))*sgn(n) 927 ENDDO 928 929 ELSE 930 ireq_mpi=ireq_mpi+1 931 buffer_r2=>message%buffers(ireq)%r2 932 !$OMP MASTER 933 CALL MPI_IRECV(buffer_r2,recv%size,MPI_REAL8,recv%rank,recv%domain+100*message%number,comm_icosa, message%mpi_req(ireq_mpi),ierr) 934 !$OMP END MASTER 935 ENDIF 936 ENDDO 937 938 ENDDO 939 940 ELSE IF (field(1)%ndim==3) THEN 941 942 ireq=0 943 ireq_mpi=0 944 DO ind=1,ndomain 945 dim3=size(field(ind)%rval3d,2) 946 rval3d=>field(ind)%rval3d 947 req=>message%request(ind) 948 949 DO isend=1,req%nsend 950 ireq=ireq+1 951 send=>req%send(isend) 952 value=>send%value 953 954 IF (send%rank/=mpi_rank .OR. .TRUE.) THEN 955 ireq_mpi=ireq_mpi+1 956 buffer_r3=>message%buffers(ireq)%r3 957 958 CALL trace_in 959 960 !$OMP DO SCHEDULE(STATIC) 961 DO n=1,send%size 962 buffer_r3(n,:)=rval3d(value(n),:) 963 ENDDO 964 965 CALL trace_out 966 967 !$OMP MASTER 968 CALL MPI_ISSEND(buffer_r3,send%size*dim3,MPI_REAL8,send%rank,ind+100*message%number,comm_icosa, message%mpi_req(ireq_mpi),ierr) 969 !$OMP END MASTER 970 ENDIF 971 ENDDO 972 973 DO irecv=1,req%nrecv 974 ireq=ireq+1 975 recv=>req%recv(irecv) 976 977 IF (recv%rank==mpi_rank .AND. .FALSE.) THEN 978 value=>recv%value 979 src_value => recv%src_value 980 src_rval3d=>field(recv%domain)%rval3d 981 sgn=>recv%sign 982 !$OMP DO SCHEDULE(STATIC) 983 DO n=1,recv%size 984 rval3d(value(n),:)=src_rval3d(src_value(n),:)*sgn(n) 985 ENDDO 986 987 ELSE 988 ireq_mpi=ireq_mpi+1 989 buffer_r3=>message%buffers(ireq)%r3 990 !$OMP MASTER 991 CALL MPI_IRECV(buffer_r3,recv%size*dim3,MPI_REAL8,recv%rank,recv%domain+100*message%number,comm_icosa, message%mpi_req(ireq_mpi),ierr) 992 !$OMP END MASTER 993 ENDIF 994 ENDDO 995 996 ENDDO 997 998 ELSE IF (field(1)%ndim==4) THEN 999 1000 ireq=0 1001 ireq_mpi=0 1002 DO ind=1,ndomain 1003 dim3=size(field(ind)%rval4d,2) 1004 dim4=size(field(ind)%rval4d,3) 1005 rval4d=>field(ind)%rval4d 1006 req=>message%request(ind) 1007 1008 DO isend=1,req%nsend 1009 ireq=ireq+1 1010 send=>req%send(isend) 1011 value=>send%value 1012 1013 IF (send%rank/=mpi_rank .OR. .TRUE.) THEN 1014 ireq_mpi=ireq_mpi+1 1015 buffer_r4=>message%buffers(ireq)%r4 1016 CALL trace_in 1017 1018 !$OMP DO SCHEDULE(STATIC) 1019 DO n=1,send%size 1020 buffer_r4(n,:,:)=rval4d(value(n),:,:) 1021 ENDDO 1022 1023 CALL trace_out 1024 1025 !$OMP MASTER 1026 CALL MPI_ISSEND(buffer_r4,send%size*dim3*dim4,MPI_REAL8,send%rank,ind+100*message%number,comm_icosa, message%mpi_req(ireq_mpi),ierr) 1027 !$OMP END MASTER 1028 ENDIF 1029 ENDDO 1030 1031 DO irecv=1,req%nrecv 1032 ireq=ireq+1 1033 recv=>req%recv(irecv) 1034 IF (recv%rank==mpi_rank .AND. .FALSE.) THEN 1035 value=>recv%value 1036 src_value => recv%src_value 1037 src_rval4d=>field(recv%domain)%rval4d 1038 sgn=>recv%sign 1039 1040 !$OMP DO SCHEDULE(STATIC) 1041 DO n=1,recv%size 1042 rval4d(value(n),:,:)=src_rval4d(src_value(n),:,:)*sgn(n) 1043 ENDDO 1044 1045 ELSE 1046 ireq_mpi=ireq_mpi+1 1047 buffer_r4=>message%buffers(ireq)%r4 1048 !$OMP MASTER 1049 CALL MPI_IRECV(buffer_r4,recv%size*dim3*dim4,MPI_REAL8,recv%rank,recv%domain+100*message%number,comm_icosa, message%mpi_req(ireq_mpi),ierr) 1050 !$OMP END MASTER 1051 ENDIF 1052 ENDDO 1053 1054 ENDDO 1055 1056 ENDIF 1057 1058 ENDIF 1059 IF (ireq_mpi /= message%nreq ) THEN 1060 STOP "ireq_mpi /= message%nreq" 1061 ENDIF 1062 1063 CALL trace_end("transfert_mpi") 1064 !$OMP BARRIER 1065 1066 END SUBROUTINE send_message_mpi 1067 1068 SUBROUTINE test_message_mpi(message) 1069 IMPLICIT NONE 1070 TYPE(t_message) :: message 1071 1072 INTEGER :: ierr 1073 !$OMP MASTER 1074 IF (.NOT. message%pending) RETURN 1075 !$OMP END MASTER 1076 1077 !$OMP MASTER 1078 IF (.NOT. message%completed) CALL MPI_TESTALL(message%nreq,message%mpi_req,message%completed,message%status,ierr) 1079 !$OMP END MASTER 1080 END SUBROUTINE test_message_mpi 1081 1082 1083 SUBROUTINE wait_message_mpi(message) 1084 USE field_mod 1085 USE domain_mod 1086 USE mpi_mod 1087 USE mpipara 1088 USE omp_para 1089 USE trace 1090 IMPLICIT NONE 1091 TYPE(t_message) :: message 1092 1093 TYPE(t_field),POINTER :: field(:) 1094 REAL(rstd),POINTER :: rval2d(:) 1095 REAL(rstd),POINTER :: rval3d(:,:) 1096 REAL(rstd),POINTER :: rval4d(:,:,:) 1097 REAL(rstd),POINTER :: buffer_r2(:) 1098 REAL(rstd),POINTER :: buffer_r3(:,:) 1099 REAL(rstd),POINTER :: buffer_r4(:,:,:) 1074 REAL(rstd),POINTER :: buffer_r(:) 1100 1075 INTEGER,POINTER :: value(:) 1101 1076 INTEGER,POINTER :: sgn(:) … … 1106 1081 INTEGER :: irecv,isend 1107 1082 INTEGER :: ireq,nreq 1108 INTEGER :: ind,n,l,m 1109 INTEGER :: dim3,dim4 1083 INTEGER :: ind,i,n,l,m 1084 INTEGER :: dim3,dim4,d3,d4 1085 INTEGER,POINTER :: src_value(:) 1086 INTEGER,POINTER :: sign(:) 1087 INTEGER :: offset,msize,rank 1088 1089 CALL trace_start("transfert_mpi") 1110 1090 1111 1091 !$OMP BARRIER 1112 1092 1113 CALL trace_start("transfert_mpi") 1114 1115 IF (.NOT. message%pending) RETURN 1116 1117 field=>message%field 1118 nreq=message%nreq 1119 1093 1094 !$OMP MASTER 1095 message%field=>field 1096 1097 IF (message%nreq>0) THEN 1098 message%completed=.FALSE. 1099 message%pending=.TRUE. 1100 ELSE 1101 message%completed=.TRUE. 1102 message%pending=.FALSE. 1103 ENDIF 1104 !$OMP END MASTER 1105 !$OMP BARRIER 1106 1120 1107 IF (field(1)%data_type==type_real) THEN 1121 1108 IF (field(1)%ndim==2) THEN 1122 1109 1123 !$OMP MASTER1124 IF (.NOT. message%completed) CALL MPI_WAITALL(nreq,message%mpi_req,message%status,ierr)1125 !$OMP END MASTER1126 !$OMP BARRIER1127 1128 ireq=01129 1110 DO ind=1,ndomain 1111 IF (.NOT. assigned_domain(ind)) CYCLE 1112 1130 1113 rval2d=>field(ind)%rval2d 1114 1131 1115 req=>message%request(ind) 1132 1133 1116 DO isend=1,req%nsend 1134 ireq=ireq+1 1135 ENDDO 1136 1137 DO irecv=1,req%nrecv 1138 ireq=ireq+1 1139 recv=>req%recv(irecv) 1140 IF (recv%rank/=mpi_rank .OR. .TRUE.) THEN 1141 buffer_r2=>message%buffers(ireq)%r2 1142 value=>recv%value 1143 sgn=>recv%sign 1144 1145 CALL trace_in 1117 send=>req%send(isend) 1118 value=>send%value 1119 1146 1120 1147 !$OMP DO SCHEDULE(STATIC) 1148 DO n=1,recv%size 1149 rval2d(value(n))=buffer_r2(n)*sgn(n) 1150 ENDDO 1151 1152 CALL trace_out 1121 IF (send%rank/=mpi_rank) THEN 1122 ireq=send%ireq 1123 1124 buffer_r=>message%buffers(ireq)%r 1125 offset=send%offset 1126 1127 DO n=1,send%size 1128 buffer_r(offset+n)=rval2d(value(n)) 1129 ENDDO 1130 1131 IF (mpi_threading_mode==MPI_THREAD_SERIALIZED) THEN 1132 !$OMP CRITICAL 1133 CALL MPI_ISEND(buffer_r,send%size,MPI_REAL8,send%rank,send%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr) 1134 !$OMP END CRITICAL 1135 ELSE IF (mpi_threading_mode==MPI_THREAD_MULTIPLE) THEN 1136 CALL MPI_ISEND(buffer_r,send%size,MPI_REAL8,send%rank,send%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr) 1137 ENDIF 1138 1153 1139 ENDIF 1154 1140 ENDDO 1155 1156 ENDDO 1157 1141 ENDDO 1142 1143 DO ind=1,ndomain 1144 IF (.NOT. assigned_domain(ind)) CYCLE 1145 rval2d=>field(ind)%rval2d 1146 req=>message%request(ind) 1147 1148 DO irecv=1,req%nrecv 1149 recv=>req%recv(irecv) 1150 1151 IF (recv%rank==mpi_rank) THEN 1152 1153 value=>recv%value 1154 src_value => recv%src_value 1155 src_rval2d=>field(recv%domain)%rval2d 1156 sgn=>recv%sign 1157 1158 DO n=1,recv%size 1159 rval2d(value(n))=src_rval2d(src_value(n))*sgn(n) 1160 ENDDO 1161 1162 1163 ELSE 1164 1165 ireq=recv%ireq 1166 buffer_r=>message%buffers(ireq)%r 1167 IF (mpi_threading_mode==MPI_THREAD_SERIALIZED) THEN 1168 !$OMP CRITICAL 1169 CALL MPI_IRECV(buffer_r,recv%size,MPI_REAL8,recv%rank,recv%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr) 1170 !$OMP END CRITICAL 1171 ELSE IF (mpi_threading_mode==MPI_THREAD_MULTIPLE) THEN 1172 CALL MPI_IRECV(buffer_r,recv%size,MPI_REAL8,recv%rank,recv%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr) 1173 ENDIF 1174 1175 ENDIF 1176 ENDDO 1177 1178 ENDDO 1179 1158 1180 1159 1181 ELSE IF (field(1)%ndim==3) THEN 1160 1161 !$OMP MASTER 1162 IF (.NOT. message%completed) CALL MPI_WAITALL(nreq,message%mpi_req,message%status,ierr) 1163 !$OMP END MASTER 1164 !$OMP BARRIER 1165 1166 ireq=0 1182 1167 1183 DO ind=1,ndomain 1184 IF (.NOT. assigned_domain(ind)) CYCLE 1185 1186 dim3=size(field(ind)%rval3d,2) 1168 1187 rval3d=>field(ind)%rval3d 1169 1188 req=>message%request(ind) 1170 1189 1171 1190 DO isend=1,req%nsend 1172 ireq=ireq+1 1173 ENDDO 1174 1175 DO irecv=1,req%nrecv 1176 ireq=ireq+1 1177 recv=>req%recv(irecv) 1178 IF (recv%rank/=mpi_rank .OR. .TRUE.) THEN 1179 buffer_r3=>message%buffers(ireq)%r3 1180 value=>recv%value 1181 sgn=>recv%sign 1182 1183 CALL trace_in 1184 1185 !$OMP DO SCHEDULE(STATIC) 1186 DO n=1,recv%size 1187 rval3d(value(n),:)=buffer_r3(n,:)*sgn(n) 1188 ENDDO 1189 1190 CALL trace_out 1191 send=>req%send(isend) 1192 value=>send%value 1193 1194 IF (send%rank/=mpi_rank) THEN 1195 ireq=send%ireq 1196 buffer_r=>message%buffers(ireq)%r 1197 offset=send%offset*dim3 1198 1199 DO d3=1,dim3 1200 DO n=1,send%size 1201 buffer_r(n+offset)=rval3d(value(n),d3) 1202 ENDDO 1203 offset=offset+send%size 1204 ENDDO 1205 1206 IF (mpi_threading_mode==MPI_THREAD_SERIALIZED) THEN 1207 !$OMP CRITICAL 1208 CALL MPI_ISEND(buffer_r,send%size*dim3,MPI_REAL8,send%rank,send%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr) 1209 !$OMP END CRITICAL 1210 ELSE IF (mpi_threading_mode==MPI_THREAD_MULTIPLE) THEN 1211 CALL MPI_ISEND(buffer_r,send%size*dim3,MPI_REAL8,send%rank,send%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr) 1212 ENDIF 1191 1213 ENDIF 1192 1214 ENDDO 1215 ENDDO 1216 1217 DO ind=1,ndomain 1218 IF (.NOT. assigned_domain(ind)) CYCLE 1219 dim3=size(field(ind)%rval3d,2) 1220 rval3d=>field(ind)%rval3d 1221 req=>message%request(ind) 1222 1223 DO irecv=1,req%nrecv 1224 recv=>req%recv(irecv) 1225 1226 IF (recv%rank==mpi_rank) THEN 1227 value=>recv%value 1228 src_value => recv%src_value 1229 src_rval3d=>field(recv%domain)%rval3d 1230 sgn=>recv%sign 1231 1232 DO n=1,recv%size 1233 rval3d(value(n),:)=src_rval3d(src_value(n),:)*sgn(n) 1234 ENDDO 1235 1236 ELSE 1237 ireq=recv%ireq 1238 buffer_r=>message%buffers(ireq)%r 1239 1240 IF (mpi_threading_mode==MPI_THREAD_SERIALIZED) THEN 1241 !$OMP CRITICAL 1242 CALL MPI_IRECV(buffer_r,recv%size*dim3,MPI_REAL8,recv%rank,recv%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr) 1243 !$OMP END CRITICAL 1244 ELSE IF (mpi_threading_mode==MPI_THREAD_MULTIPLE) THEN 1245 CALL MPI_IRECV(buffer_r,recv%size*dim3,MPI_REAL8,recv%rank,recv%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr) 1246 ENDIF 1247 ENDIF 1248 ENDDO 1193 1249 1194 1250 ENDDO 1195 1251 1196 1252 ELSE IF (field(1)%ndim==4) THEN 1197 !$OMP MASTER 1198 IF (.NOT. message%completed) CALL MPI_WAITALL(nreq,message%mpi_req,message%status,ierr) 1199 !$OMP END MASTER 1200 !$OMP BARRIER 1201 1202 ireq=0 1253 1203 1254 DO ind=1,ndomain 1255 IF (.NOT. assigned_domain(ind)) CYCLE 1256 1257 dim3=size(field(ind)%rval4d,2) 1258 dim4=size(field(ind)%rval4d,3) 1204 1259 rval4d=>field(ind)%rval4d 1205 1260 req=>message%request(ind) 1206 1261 1207 1262 DO isend=1,req%nsend 1208 ireq=ireq+11209 ENDDO1210 1211 DO irecv=1,req%nrecv1212 ireq=ireq+1 1213 recv=>req%recv(irecv)1214 IF (recv%rank/=mpi_rank .OR. .TRUE.) THEN1215 buffer_r4=>message%buffers(ireq)%r41216 value=>recv%value 1217 sgn=>recv%sign1218 1219 CALL trace_in1220 1221 !$OMP DO SCHEDULE(STATIC) 1222 DO n=1,recv%size1223 rval4d(value(n),:,:)=buffer_r4(n,:,:)*sgn(n)1263 send=>req%send(isend) 1264 value=>send%value 1265 1266 IF (send%rank/=mpi_rank) THEN 1267 1268 ireq=send%ireq 1269 buffer_r=>message%buffers(ireq)%r 1270 offset=send%offset*dim3*dim4 1271 1272 DO d4=1,dim4 1273 DO d3=1,dim3 1274 DO n=1,send%size 1275 buffer_r(n+offset)=rval4d(value(n),d3,d4) 1276 ENDDO 1277 offset=offset+send%size 1278 ENDDO 1224 1279 ENDDO 1225 1280 1226 CALL trace_out 1281 IF (mpi_threading_mode==MPI_THREAD_SERIALIZED) THEN 1282 !$OMP CRITICAL 1283 CALL MPI_ISEND(buffer_r,send%size*dim3*dim4,MPI_REAL8,send%rank,send%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr) 1284 !$OMP END CRITICAL 1285 ELSE IF (mpi_threading_mode==MPI_THREAD_MULTIPLE) THEN 1286 CALL MPI_ISEND(buffer_r,send%size*dim3*dim4,MPI_REAL8,send%rank,send%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr) 1287 ENDIF 1288 1227 1289 ENDIF 1228 1290 ENDDO 1229 1230 ENDDO 1231 1291 ENDDO 1292 1293 DO ind=1,ndomain 1294 IF (.NOT. assigned_domain(ind)) CYCLE 1295 1296 dim3=size(field(ind)%rval4d,2) 1297 dim4=size(field(ind)%rval4d,3) 1298 rval4d=>field(ind)%rval4d 1299 req=>message%request(ind) 1300 1301 DO irecv=1,req%nrecv 1302 recv=>req%recv(irecv) 1303 IF (recv%rank==mpi_rank) THEN 1304 value=>recv%value 1305 src_value => recv%src_value 1306 src_rval4d=>field(recv%domain)%rval4d 1307 sgn=>recv%sign 1308 1309 DO n=1,recv%size 1310 rval4d(value(n),:,:)=src_rval4d(src_value(n),:,:)*sgn(n) 1311 ENDDO 1312 1313 ELSE 1314 1315 ireq=recv%ireq 1316 buffer_r=>message%buffers(ireq)%r 1317 IF (mpi_threading_mode==MPI_THREAD_SERIALIZED) THEN 1318 !$OMP CRITICAL 1319 CALL MPI_IRECV(buffer_r,recv%size*dim3*dim4,MPI_REAL8,recv%rank,recv%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr) 1320 !$OMP END CRITICAL 1321 ELSE IF (mpi_threading_mode==MPI_THREAD_MULTIPLE) THEN 1322 CALL MPI_IRECV(buffer_r,recv%size*dim3*dim4,MPI_REAL8,recv%rank,recv%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr) 1323 ENDIF 1324 1325 ENDIF 1326 ENDDO 1327 ENDDO 1328 1232 1329 ENDIF 1233 1330 1331 IF (mpi_threading_mode==MPI_THREAD_FUNNELED) THEN 1332 !$OMP BARRIER 1333 !$OMP MASTER 1334 1335 DO ireq=1,message%nreq_send 1336 buffer_r=>message%buffers(ireq)%r 1337 msize=message%buffers(ireq)%size 1338 rank=message%buffers(ireq)%rank 1339 CALL MPI_ISEND(buffer_r,msize,MPI_REAL8,rank,1000*message%number,comm_icosa, message%mpi_req(ireq),ierr) 1340 ENDDO 1341 1342 DO ireq=message%nreq_send+1,message%nreq_send+message%nreq_recv 1343 buffer_r=>message%buffers(ireq)%r 1344 msize=message%buffers(ireq)%size 1345 rank=message%buffers(ireq)%rank 1346 CALL MPI_IRECV(buffer_r,msize,MPI_REAL8,rank,1000*message%number,comm_icosa, message%mpi_req(ireq),ierr) 1347 ENDDO 1348 1349 !$OMP END MASTER 1350 ENDIF 1234 1351 ENDIF 1352 1353 !$OMP BARRIER 1354 CALL trace_end("transfert_mpi") 1355 1356 END SUBROUTINE send_message_mpi 1357 1358 SUBROUTINE test_message_mpi(message) 1359 IMPLICIT NONE 1360 TYPE(t_message) :: message 1361 1362 INTEGER :: ierr 1235 1363 1236 1364 !$OMP MASTER 1237 message%pending=.FALSE.1365 IF (message%pending .AND. .NOT. message%completed) CALL MPI_TESTALL(message%nreq,message%mpi_req,message%completed,message%status,ierr) 1238 1366 !$OMP END MASTER 1239 1240 CALL trace_end("transfert_mpi") 1241 !$OMP BARRIER 1242 1243 END SUBROUTINE wait_message_mpi 1244 1245 1246 SUBROUTINE transfert_request_mpi(field,request) 1367 END SUBROUTINE test_message_mpi 1368 1369 1370 SUBROUTINE wait_message_mpi(message) 1247 1371 USE field_mod 1248 1372 USE domain_mod 1249 1373 USE mpi_mod 1250 1374 USE mpipara 1375 USE omp_para 1251 1376 USE trace 1252 1377 IMPLICIT NONE 1378 TYPE(t_message) :: message 1379 1253 1380 TYPE(t_field),POINTER :: field(:) 1254 TYPE(t_request),POINTER :: request(:)1255 1381 REAL(rstd),POINTER :: rval2d(:) 1256 1382 REAL(rstd),POINTER :: rval3d(:,:) 1257 1383 REAL(rstd),POINTER :: rval4d(:,:,:) 1258 REAL(rstd),POINTER :: buffer_r2(:) 1259 REAL(rstd),POINTER :: buffer_r3(:,:) 1260 REAL(rstd),POINTER :: buffer_r4(:,:,:) 1384 REAL(rstd),POINTER :: buffer_r(:) 1261 1385 INTEGER,POINTER :: value(:) 1262 1386 INTEGER,POINTER :: sgn(:) … … 1267 1391 INTEGER :: irecv,isend 1268 1392 INTEGER :: ireq,nreq 1269 INTEGER :: ind,n 1270 INTEGER :: dim3,dim4 1393 INTEGER :: ind,n,l,m,i 1394 INTEGER :: dim3,dim4,d3,d4 1395 INTEGER :: offset 1396 1397 IF (.NOT. message%pending) RETURN 1271 1398 1272 1399 CALL trace_start("transfert_mpi") 1273 1400 1401 field=>message%field 1402 nreq=message%nreq 1403 1274 1404 IF (field(1)%data_type==type_real) THEN 1275 1405 IF (field(1)%ndim==2) THEN 1276 1277 nreq=sum(request(:)%nsend)+sum(request(:)%nrecv) 1278 ALLOCATE(mpi_req(nreq))1279 ALLOCATE(status(MPI_STATUS_SIZE,nreq)) 1280 1281 ireq=01406 1407 !$OMP MASTER 1408 IF (.NOT. message%completed) CALL MPI_WAITALL(nreq,message%mpi_req,message%status,ierr) 1409 !$OMP END MASTER 1410 !$OMP BARRIER 1411 1282 1412 DO ind=1,ndomain 1413 IF (.NOT. assigned_domain(ind)) CYCLE 1414 1283 1415 rval2d=>field(ind)%rval2d 1284 1285 req=>request(ind) 1286 DO isend=1,req%nsend 1287 send=>req%send(isend) 1288 1289 ALLOCATE(send%buffer_r2(send%size)) 1290 buffer_r2=>send%buffer_r2 1291 value=>send%value 1292 DO n=1,send%size 1293 buffer_r2(n)=rval2d(value(n)) 1294 ENDDO 1295 1296 ireq=ireq+1 1297 CALL MPI_ISEND(buffer_r2,send%size,MPI_REAL8,send%rank,ind,comm_icosa, mpi_req(ireq),ierr) 1298 ENDDO 1299 1416 req=>message%request(ind) 1300 1417 DO irecv=1,req%nrecv 1301 1418 recv=>req%recv(irecv) 1302 ALLOCATE(recv%buffer_r2(recv%size)) 1303 1304 ireq=ireq+1 1305 CALL MPI_IRECV(recv%buffer_r2,recv%size,MPI_REAL8,recv%rank,recv%domain,comm_icosa, mpi_req(ireq),ierr) 1419 IF (recv%rank/=mpi_rank) THEN 1420 ireq=recv%ireq 1421 buffer_r=>message%buffers(ireq)%r 1422 value=>recv%value 1423 sgn=>recv%sign 1424 offset=recv%offset 1425 DO n=1,recv%size 1426 rval2d(value(n))=buffer_r(n+offset)*sgn(n) 1427 ENDDO 1428 1429 ENDIF 1306 1430 ENDDO 1307 1431 1308 1432 ENDDO 1309 1310 CALL MPI_WAITALL(nreq,mpi_req,status,ierr) 1433 1434 1435 ELSE IF (field(1)%ndim==3) THEN 1436 1437 !$OMP MASTER 1438 IF (.NOT. message%completed) CALL MPI_WAITALL(nreq,message%mpi_req,message%status,ierr) 1439 !$OMP END MASTER 1440 !$OMP BARRIER 1441 1311 1442 1312 1443 DO ind=1,ndomain 1313 rval2d=>field(ind)%rval2d 1314 1315 req=>request(ind) 1316 DO isend=1,req%nsend 1317 send=>req%send(isend) 1318 DEALLOCATE(send%buffer_r2) 1319 ENDDO 1320 1444 IF (.NOT. assigned_domain(ind)) CYCLE 1445 1446 rval3d=>field(ind)%rval3d 1447 req=>message%request(ind) 1321 1448 DO irecv=1,req%nrecv 1322 1449 recv=>req%recv(irecv) 1323 buffer_r2=>recv%buffer_r2 1324 value=>recv%value 1325 sgn=>recv%sign 1326 DO n=1,recv%size 1327 rval2d(value(n))=buffer_r2(n)*sgn(n) 1328 ENDDO 1329 DEALLOCATE(recv%buffer_r2) 1450 IF (recv%rank/=mpi_rank) THEN 1451 ireq=recv%ireq 1452 buffer_r=>message%buffers(ireq)%r 1453 value=>recv%value 1454 sgn=>recv%sign 1455 1456 dim3=size(rval3d,2) 1457 offset=recv%offset*dim3 1458 DO d3=1,dim3 1459 DO n=1,recv%size 1460 rval3d(value(n),d3)=buffer_r(n+offset)*sgn(n) 1461 ENDDO 1462 offset=offset+recv%size 1463 ENDDO 1464 ENDIF 1330 1465 ENDDO 1331 1466 1332 1467 ENDDO 1333 1334 1335 ELSE IF (field(1)%ndim==3) THEN 1336 1337 nreq=sum(request(:)%nsend)+sum(request(:)%nrecv) 1338 ALLOCATE(mpi_req(nreq)) 1339 ALLOCATE(status(MPI_STATUS_SIZE,nreq)) 1340 1341 ireq=0 1468 1469 ELSE IF (field(1)%ndim==4) THEN 1470 !$OMP MASTER 1471 IF (.NOT. message%completed) CALL MPI_WAITALL(nreq,message%mpi_req,message%status,ierr) 1472 !$OMP END MASTER 1473 !$OMP BARRIER 1474 1475 1342 1476 DO ind=1,ndomain 1343 dim3=size(field(ind)%rval3d,2) 1344 rval3d=>field(ind)%rval3d 1345 1346 req=>request(ind) 1347 DO isend=1,req%nsend 1348 send=>req%send(isend) 1349 1350 ALLOCATE(send%buffer_r3(send%size,dim3)) 1351 buffer_r3=>send%buffer_r3 1352 value=>send%value 1353 DO n=1,send%size 1354 buffer_r3(n,:)=rval3d(value(n),:) 1355 ENDDO 1356 1357 ireq=ireq+1 1358 CALL MPI_ISEND(buffer_r3,send%size*dim3,MPI_REAL8,send%rank,ind,comm_icosa, mpi_req(ireq),ierr) 1359 ENDDO 1360 1477 IF (.NOT. assigned_domain(ind)) CYCLE 1478 1479 rval4d=>field(ind)%rval4d 1480 req=>message%request(ind) 1361 1481 DO irecv=1,req%nrecv 1362 1482 recv=>req%recv(irecv) 1363 ALLOCATE(recv%buffer_r3(recv%size,dim3)) 1364 1365 ireq=ireq+1 1366 CALL MPI_IRECV(recv%buffer_r3,recv%size*dim3,MPI_REAL8,recv%rank,recv%domain,comm_icosa, mpi_req(ireq),ierr) 1483 IF (recv%rank/=mpi_rank) THEN 1484 ireq=recv%ireq 1485 buffer_r=>message%buffers(ireq)%r 1486 value=>recv%value 1487 sgn=>recv%sign 1488 1489 dim3=size(rval4d,2) 1490 dim4=size(rval4d,3) 1491 offset=recv%offset*dim3*dim4 1492 DO d4=1,dim4 1493 DO d3=1,dim3 1494 DO n=1,recv%size 1495 rval4d(value(n),d3,d4)=buffer_r(n+offset)*sgn(n) 1496 ENDDO 1497 offset=offset+recv%size 1498 ENDDO 1499 ENDDO 1500 ENDIF 1367 1501 ENDDO 1368 1502 1369 1503 ENDDO 1370 1371 CALL MPI_WAITALL(nreq,mpi_req,status,ierr)1372 1373 DO ind=1,ndomain1374 rval3d=>field(ind)%rval3d1375 1376 req=>request(ind)1377 DO isend=1,req%nsend1378 send=>req%send(isend)1379 DEALLOCATE(send%buffer_r3)1380 ENDDO1381 1382 DO irecv=1,req%nrecv1383 recv=>req%recv(irecv)1384 buffer_r3=>recv%buffer_r31385 value=>recv%value1386 sgn=>recv%sign1387 DO n=1,recv%size1388 rval3d(value(n),:)=buffer_r3(n,:)*sgn(n)1389 ENDDO1390 DEALLOCATE(recv%buffer_r3)1391 ENDDO1392 1393 ENDDO1394 1395 ELSE IF (field(1)%ndim==4) THEN1396 1397 nreq=sum(request(:)%nsend)+sum(request(:)%nrecv)1398 ALLOCATE(mpi_req(nreq))1399 ALLOCATE(status(MPI_STATUS_SIZE,nreq))1400 1401 ireq=01402 DO ind=1,ndomain1403 dim3=size(field(ind)%rval4d,2)1404 dim4=size(field(ind)%rval4d,3)1405 rval4d=>field(ind)%rval4d1406 1407 req=>request(ind)1408 DO isend=1,req%nsend1409 send=>req%send(isend)1410 1411 ALLOCATE(send%buffer_r4(send%size,dim3,dim4))1412 buffer_r4=>send%buffer_r41413 value=>send%value1414 DO n=1,send%size1415 buffer_r4(n,:,:)=rval4d(value(n),:,:)1416 ENDDO1417 1418 ireq=ireq+11419 CALL MPI_ISEND(buffer_r4,send%size*dim3*dim4,MPI_REAL8,send%rank,ind,comm_icosa, mpi_req(ireq),ierr)1420 ENDDO1421 1422 DO irecv=1,req%nrecv1423 recv=>req%recv(irecv)1424 ALLOCATE(recv%buffer_r4(recv%size,dim3,dim4))1425 1426 ireq=ireq+11427 CALL MPI_IRECV(recv%buffer_r4,recv%size*dim3*dim4,MPI_REAL8,recv%rank,recv%domain,comm_icosa, mpi_req(ireq),ierr)1428 ENDDO1429 1430 ENDDO1431 1432 CALL MPI_WAITALL(nreq,mpi_req,status,ierr)1433 1434 DO ind=1,ndomain1435 rval4d=>field(ind)%rval4d1436 1437 req=>request(ind)1438 DO isend=1,req%nsend1439 send=>req%send(isend)1440 DEALLOCATE(send%buffer_r4)1441 ENDDO1442 1443 DO irecv=1,req%nrecv1444 recv=>req%recv(irecv)1445 buffer_r4=>recv%buffer_r41446 value=>recv%value1447 sgn=>recv%sign1448 DO n=1,recv%size1449 rval4d(value(n),:,:)=buffer_r4(n,:,:)*sgn(n)1450 ENDDO1451 DEALLOCATE(recv%buffer_r4)1452 ENDDO1453 1454 ENDDO1455 1504 1456 1505 ENDIF … … 1458 1507 ENDIF 1459 1508 1509 !$OMP MASTER 1510 message%pending=.FALSE. 1511 !$OMP END MASTER 1512 1460 1513 CALL trace_end("transfert_mpi") 1461 1514 !$OMP BARRIER 1515 1516 END SUBROUTINE wait_message_mpi 1517 1518 SUBROUTINE transfert_request_mpi(field,request) 1519 USE field_mod 1520 IMPLICIT NONE 1521 TYPE(t_field),POINTER :: field(:) 1522 TYPE(t_request),POINTER :: request(:) 1523 1524 TYPE(t_message),SAVE :: message 1525 1526 1527 CALL init_message_mpi(field,request, message) 1528 CALL transfert_message_mpi(field,message) 1529 CALL finalize_message_mpi(field,message) 1530 1462 1531 END SUBROUTINE transfert_request_mpi 1532 1533 1463 1534 1464 1535 SUBROUTINE transfert_request_seq(field,request) -
codes/icosagcm/trunk/src/vertical_interp.f90
r105 r186 3 3 PRIVATE 4 4 5 TYPE(t_field), POINTER :: f_p(:)5 TYPE(t_field),SAVE, POINTER :: f_p(:) 6 6 7 7 … … 36 36 37 37 DO ind=1,ndomain 38 IF (.NOT. assigned_domain(ind)) CYCLE 38 39 CALL swap_dimensions(ind) 39 40 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/vorticity.f90
r146 r186 16 16 17 17 DO ind=1,ndomain 18 IF (.NOT. assigned_domain(ind)) CYCLE 18 19 CALL swap_dimensions(ind) 19 20 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/wind.f90
r179 r186 13 13 14 14 DO ind=1,ndomain 15 IF (.NOT. assigned_domain(ind)) CYCLE 15 16 CALL swap_dimensions(ind) 16 17 CALL swap_geometry(ind) -
codes/icosagcm/trunk/src/write_field.f90
r171 r186 1770 1770 IMPLICIT NONE 1771 1771 INTEGER :: i,k,status 1772 1772 !$OMP MASTER 1773 1773 DO i=1,NbField 1774 1774 status=NF90_CLOSE(FieldId(i)) 1775 1775 ENDDO 1776 !$OMP END MASTER 1776 1777 END SUBROUTINE Close_files 1777 1778 -
codes/icosagcm/trunk/src/xios_mod.F90
r173 r186 7 7 PUBLIC 8 8 LOGICAL,SAVE :: using_xios 9 9 10 INTEGER,SAVE :: ncell_i 11 !$OMP THREADPRIVATE(ncell_i) 10 12 INTEGER,SAVE :: ncell_v 13 !$OMP THREADPRIVATE(ncell_v) 11 14 12 15 PRIVATE ncell_i,ncell_v … … 17 20 18 21 SUBROUTINE xios_init 19 USE IOIPSL22 USE getin_mod 20 23 IMPLICIT NONE 21 24 … … 46 49 TYPE(t_domain),POINTER :: d 47 50 51 !$OMP BARRIER 52 !$OMP MASTER 48 53 CALL xios_context_initialize("icosagcm",comm_icosa) 49 54 CALL xios_get_handle("icosagcm",ctx_hdl) … … 186 191 CALL xios_set_timestep(dtime) 187 192 CALL xios_close_context_definition() 188 193 !$OMP END MASTER 194 !$OMP BARRIER 195 189 196 END SUBROUTINE xios_init_write_field 190 197 … … 197 204 CHARACTER(LEN=10) :: str_number 198 205 INTEGER :: iq 206 207 !$OMP BARRIER 208 !$OMP MASTER 199 209 200 210 IF (Field(1)%field_type==field_T) THEN … … 225 235 ENDIF 226 236 ENDIF 237 !$OMP END MASTER 238 !$OMP BARRIER 227 239 228 240 END SUBROUTINE xios_write_field … … 400 412 SUBROUTINE xios_write_field_finalize 401 413 IMPLICIT NONE 414 415 !$OMP BARRIER 416 !$OMP MASTER 402 417 CALL xios_context_finalize 418 !$OMP END MASTER 419 !$OMP BARRIER 420 403 421 END SUBROUTINE xios_write_field_finalize 404 422
Note: See TracChangeset
for help on using the changeset viewer.