Changeset 845 for codes/icosagcm/devel/src
- Timestamp:
- 05/04/19 21:35:04 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/dissip/dissip_gcm.f90
r714 r845 1 1 MODULE dissip_gcm_mod 2 2 USE icosa 3 3 USE omp_para 4 USE trace 5 IMPLICIT NONE 4 6 PRIVATE 5 7 … … 41 43 REAL,SAVE :: dtdissip 42 44 !$OMP THREADPRIVATE(dtdissip) 43 45 44 46 PUBLIC init_dissip, dissip 47 45 48 CONTAINS 46 49 47 50 SUBROUTINE allocate_dissip 48 USE icosa49 IMPLICIT NONE50 51 CALL allocate_field(f_due_diss1,field_u,type_real,llm) 51 52 CALL allocate_field(f_due_diss2,field_u,type_real,llm) … … 58 59 59 60 SUBROUTINE init_dissip 60 USE icosa 61 USE disvert_mod 62 USE mpi_mod 63 USE mpipara 64 USE transfert_mod 65 USE time_mod 66 USE transfert_omp_mod 67 USE omp_para 68 IMPLICIT NONE 69 70 TYPE(t_field),POINTER,SAVE :: f_u(:) 71 TYPE(t_field),POINTER,SAVE :: f_du(:) 72 REAL(rstd),POINTER :: u(:) 73 REAL(rstd),POINTER :: du(:) 74 TYPE(t_field),POINTER,SAVE :: f_theta(:) 75 TYPE(t_field),POINTER ,SAVE :: f_dtheta(:) 76 REAL(rstd),POINTER :: theta(:) 77 REAL(rstd),POINTER :: dtheta(:) 78 REAL(rstd) :: dumax,dumax1 79 REAL(rstd) :: dthetamax,dthetamax1 80 REAL(rstd) :: r 81 REAL(rstd) :: tau 82 REAL(rstd) :: zz, zvert, fact 83 INTEGER :: l 84 CHARACTER(len=255) :: rayleigh_friction_key 85 REAL(rstd) :: mintau 86 INTEGER :: seed_size 87 INTEGER,ALLOCATABLE :: seed(:) 88 89 90 INTEGER :: i,j,ij,ind,it,iter,M 91 92 rayleigh_friction_key='none' 93 CALL getin("rayleigh_friction_type",rayleigh_friction_key) 94 SELECT CASE(TRIM(rayleigh_friction_key)) 95 CASE('none') 96 rayleigh_friction_type=0 97 IF (is_master) PRINT *, 'No Rayleigh friction' 98 CASE('dcmip2_schaer_noshear') 99 rayleigh_friction_type=1 100 rayleigh_shear=0 101 IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain without shear DCMIP2.1' 102 CASE('dcmip2_schaer_shear') 103 rayleigh_shear=1 104 rayleigh_friction_type=2 105 IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2' 106 CASE('giant_liu_schneider') 107 rayleigh_friction_type=99 108 IF (is_master) PRINT *, 'Rayleigh friction : giant planets Liu Schneider 2010' 109 CASE DEFAULT 110 IF (is_master) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), & 111 ' in dissip_gcm.f90/init_dissip' 112 STOP 113 END SELECT 114 115 IF(rayleigh_friction_type>0) THEN 116 rayleigh_tau=0. 117 CALL getin("rayleigh_friction_tau",rayleigh_tau) 118 rayleigh_tau = rayleigh_tau / scale_factor 119 IF(rayleigh_tau<=0) THEN 120 IF (is_master) PRINT *, 'Forbidden : negative value for rayleigh_friction_tau =',rayleigh_tau 121 STOP 122 END IF 123 IF(rayleigh_friction_type == 99) THEN 124 rayleigh_limlat=0. 125 CALL getin("rayleigh_limlat",rayleigh_limlat) 126 rayleigh_limlat = rayleigh_limlat*3.14159/180. 127 ENDIF 128 END IF 129 61 REAL(rstd) :: tau 62 CHARACTER(len=255) :: rayleigh_friction_key 63 rayleigh_friction_key='none' 64 CALL getin("rayleigh_friction_type",rayleigh_friction_key) 65 SELECT CASE(TRIM(rayleigh_friction_key)) 66 CASE('none') 67 rayleigh_friction_type=0 68 IF (is_master) PRINT *, 'No Rayleigh friction' 69 CASE('dcmip2_schaer_noshear') 70 rayleigh_friction_type=1 71 rayleigh_shear=0 72 IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain without shear DCMIP2.1' 73 CASE('dcmip2_schaer_shear') 74 rayleigh_shear=1 75 rayleigh_friction_type=2 76 IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2' 77 CASE('giant_liu_schneider') 78 rayleigh_friction_type=99 79 IF (is_master) PRINT *, 'Rayleigh friction : giant planets Liu Schneider 2010' 80 CASE DEFAULT 81 IF (is_master) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), & 82 ' in dissip_gcm.f90/init_dissip' 83 STOP 84 END SELECT 85 86 IF(rayleigh_friction_type>0) THEN 87 rayleigh_tau=0. 88 CALL getin("rayleigh_friction_tau",rayleigh_tau) 89 rayleigh_tau = rayleigh_tau / scale_factor 90 IF(rayleigh_tau<=0) THEN 91 IF (is_master) PRINT *, 'Forbidden : negative value for rayleigh_friction_tau =',rayleigh_tau 92 STOP 93 END IF 94 IF(rayleigh_friction_type == 99) THEN 95 rayleigh_limlat=0. 96 CALL getin("rayleigh_limlat",rayleigh_limlat) 97 rayleigh_limlat = rayleigh_limlat*3.14159/180. 98 ENDIF 99 END IF 100 130 101 CALL allocate_dissip 102 103 CALL init_message(f_due_diss1,req_e1_vect,req_due) 104 CALL init_message(f_dtheta_diss,req_i1,req_dtheta) 105 106 tau_graddiv(:)=5000 107 CALL getin("tau_graddiv",tau) 108 tau_graddiv(:)=tau/scale_factor 109 110 CALL getin("nitergdiv",nitergdiv) 111 112 tau_gradrot(:)=5000 113 CALL getin("tau_gradrot",tau) 114 tau_gradrot(:)=tau/scale_factor 115 116 CALL getin("nitergrot",nitergrot) 117 118 119 tau_divgrad(:)=5000 120 CALL getin("tau_divgrad",tau) 121 tau_divgrad(:)=tau/scale_factor 122 123 CALL getin("niterdivgrad",niterdivgrad) 124 125 CALL dissip_constants 126 CALL dissip_profile 127 END SUBROUTINE init_dissip 128 129 SUBROUTINE dissip_constants 130 ! The SAVE attribute is required here, otherwise allocate_field will not work with OpenMP 131 TYPE(t_field),POINTER, SAVE :: f_u(:) 132 TYPE(t_field),POINTER, SAVE :: f_du(:) 133 TYPE(t_field),POINTER, SAVE :: f_theta(:) 134 TYPE(t_field),POINTER, SAVE :: f_dtheta(:) 135 REAL(rstd),POINTER :: u(:) 136 REAL(rstd),POINTER :: du(:) 137 REAL(rstd),POINTER :: theta(:) 138 REAL(rstd),POINTER :: dtheta(:) 139 REAL(rstd) :: dumax, dthetamax 140 INTEGER :: it, iter, ind 131 141 CALL allocate_field(f_u,field_u,type_real) 132 142 CALL allocate_field(f_du,field_u,type_real) 133 143 CALL allocate_field(f_theta,field_t,type_real) 134 144 CALL allocate_field(f_dtheta,field_t,type_real) 135 136 CALL init_message(f_due_diss1,req_e1_vect,req_due) 137 CALL init_message(f_dtheta_diss,req_i1,req_dtheta) 138 139 tau_graddiv(:)=5000 140 CALL getin("tau_graddiv",tau) 141 tau_graddiv(:)=tau/scale_factor 142 143 CALL getin("nitergdiv",nitergdiv) 144 145 tau_gradrot(:)=5000 146 CALL getin("tau_gradrot",tau) 147 tau_gradrot(:)=tau/scale_factor 148 149 CALL getin("nitergrot",nitergrot) 150 151 152 tau_divgrad(:)=5000 153 CALL getin("tau_divgrad",tau) 154 tau_divgrad(:)=tau/scale_factor 155 156 CALL getin("niterdivgrad",niterdivgrad) 157 158 159 cgraddiv=-1 160 cdivgrad=-1 161 cgradrot=-1 162 163 !$OMP BARRIER 164 !$OMP MASTER 165 DO ind=1,ndomain 166 CALL swap_dimensions(ind) 167 CALL swap_geometry(ind) 168 u=f_u(ind) 169 170 ! set random seed to get reproductibility when using a different number of process 171 CALL RANDOM_SEED(size=M) 172 CALL RANDOM_SEED(put=(/(i,i=1,M)/)) 173 174 DO j=jj_begin,jj_end 175 DO i=ii_begin,ii_end 176 ij=(j-1)*iim+i 177 CALL RANDOM_NUMBER(r) 178 u(ij+u_right)=r-0.5 179 CALL RANDOM_NUMBER(r) 180 u(ij+u_lup)=r-0.5 181 CALL RANDOM_NUMBER(r) 182 u(ij+u_ldown)=r-0.5 183 ENDDO 184 ENDDO 185 ENDDO 186 !$OMP END MASTER 187 !$OMP BARRIER 188 145 146 PRINT *, '=========', SHAPE(f_u) 147 148 !--------------------------- Compute cgraddiv ---------------------------- 149 cgraddiv=-1. 150 CALL random_vel(f_u) 189 151 DO it=1,20 190 191 dumax=0 192 DO iter=1,nitergdiv 193 CALL transfert_request(f_u,req_e1_vect) 194 DO ind=1,ndomain 195 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 196 CALL swap_dimensions(ind) 197 CALL swap_geometry(ind) 198 u=f_u(ind) 199 du=f_du(ind) 200 CALL compute_gradiv(u,du,1,1) 201 u=du 202 ENDDO 203 ENDDO 204 205 CALL transfert_request(f_du,req_e1_vect) 206 207 DO ind=1,ndomain 208 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 209 CALL swap_dimensions(ind) 210 CALL swap_geometry(ind) 211 u=f_u(ind) 212 du=f_du(ind) 213 214 DO j=jj_begin,jj_end 215 DO i=ii_begin,ii_end 216 ij=(j-1)*iim+i 217 if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right))) 218 if (le(ij+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_lup))) 219 if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown))) 152 DO iter=1,nitergdiv 153 CALL transfert_request(f_u,req_e1_vect) 154 DO ind=1,ndomain 155 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 156 CALL swap_dimensions(ind) 157 CALL swap_geometry(ind) 158 u=f_u(ind) 159 du=f_du(ind) 160 CALL compute_gradiv(u,du,1,1) 220 161 ENDDO 221 ENDDO 222 ENDDO 223 224 IF (using_mpi) THEN 225 CALL reduce_max_omp(dumax,dumax1) 226 !$OMP MASTER 227 CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 228 !$OMP END MASTER 229 CALL bcast_omp(dumax) 230 ELSE 231 CALL allreduce_max_omp(dumax,dumax1) 232 dumax=dumax1 233 ENDIF 234 235 DO ind=1,ndomain 236 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 237 CALL swap_dimensions(ind) 238 CALL swap_geometry(ind) 239 u=f_u(ind) 240 du=f_du(ind) 241 u=du/dumax 242 ENDDO 243 IF (is_master) PRINT *,"gradiv : it :",it ,": dumax",dumax 244 245 ENDDO 246 IF (is_master) PRINT *,"gradiv : dumax",dumax 247 IF (is_master) PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., & 248 'effective T-cell half-edge size (km)', dumax**(-.5/nitergdiv)/1000 249 IF (is_master) PRINT *, 'Max. time step assuming c=340 m/s and Courant number=3 (ARK2.3) :', & 250 3./340.*dumax**(-.5/nitergdiv) 162 ENDDO 163 CALL transfert_request(f_du,req_e1_vect) 164 dumax = max_vel(f_du) 165 CALL rescale_field(1./dumax, f_du, f_u) 166 IF (is_master) PRINT *,"gradiv : it :",it ,": dumax",dumax 167 ENDDO 251 168 252 169 cgraddiv=dumax**(-1./nitergdiv) 253 IF (is_master) PRINT *,"cgraddiv : ",cgraddiv 254 255 !$OMP BARRIER 256 !$OMP MASTER 257 DO ind=1,ndomain 258 CALL swap_dimensions(ind) 259 CALL swap_geometry(ind) 260 u=f_u(ind) 261 262 ! set random seed to get reproductibility when using a different number of process 263 CALL RANDOM_SEED(size=M) 264 CALL RANDOM_SEED(put=(/(i,i=1,M)/)) 265 266 DO j=jj_begin,jj_end 267 DO i=ii_begin,ii_end 268 ij=(j-1)*iim+i 269 CALL RANDOM_NUMBER(r) 270 u(ij+u_right)=r-0.5 271 CALL RANDOM_NUMBER(r) 272 u(ij+u_lup)=r-0.5 273 CALL RANDOM_NUMBER(r) 274 u(ij+u_ldown)=r-0.5 275 ENDDO 276 ENDDO 277 ENDDO 278 !$OMP END MASTER 279 !$OMP BARRIER 280 281 170 IF (is_master) THEN 171 PRINT *,"gradiv : dumax",dumax 172 PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., & 173 'effective T-cell half-edge size (km)', dumax**(-.5/nitergdiv)/1000 174 PRINT *, 'Max. time step assuming c=340 m/s and Courant number=3 (ARK2.3) :', & 175 3./340.*dumax**(-.5/nitergdiv) 176 PRINT *,"cgraddiv : ",cgraddiv 177 END IF 178 179 !----------------- Compute cgradrot -------------------- 180 cgradrot=-1. 181 CALL random_vel(f_u) 282 182 DO it=1,20 283 284 dumax=0 285 DO iter=1,nitergrot 286 CALL transfert_request(f_u,req_e1_vect) 287 DO ind=1,ndomain 288 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 289 CALL swap_dimensions(ind) 290 CALL swap_geometry(ind) 291 u=f_u(ind) 292 du=f_du(ind) 293 CALL compute_gradrot(u,du,1,1) 294 u=du 295 ENDDO 296 ENDDO 297 298 CALL transfert_request(f_du,req_e1_vect) 299 300 DO ind=1,ndomain 301 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 302 CALL swap_dimensions(ind) 303 CALL swap_geometry(ind) 304 u=f_u(ind) 305 du=f_du(ind) 306 307 DO j=jj_begin,jj_end 308 DO i=ii_begin,ii_end 309 ij=(j-1)*iim+i 310 if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right))) 311 if (le(ij+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_lup))) 312 if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown))) 183 DO iter=1,nitergrot 184 CALL transfert_request(f_u,req_e1_vect) 185 DO ind=1,ndomain 186 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 187 CALL swap_dimensions(ind) 188 CALL swap_geometry(ind) 189 u=f_u(ind) 190 du=f_du(ind) 191 CALL compute_gradrot(u,du,1,1) 192 u=du 313 193 ENDDO 314 ENDDO 315 ENDDO 316 317 IF (using_mpi) THEN 318 CALL reduce_max_omp(dumax,dumax1) 319 !$OMP MASTER 320 CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 321 !$OMP END MASTER 322 CALL bcast_omp(dumax) 323 ELSE 324 CALL allreduce_max_omp(dumax,dumax1) 325 dumax=dumax1 326 ENDIF 327 328 329 DO ind=1,ndomain 330 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 331 CALL swap_dimensions(ind) 332 CALL swap_geometry(ind) 333 u=f_u(ind) 334 du=f_du(ind) 335 u=du/dumax 336 ENDDO 337 338 IF (is_master) PRINT *,"gradrot : it :",it ,": dumax",dumax 339 340 ENDDO 341 IF (is_master) PRINT *,"gradrot : dumax",dumax 342 194 ENDDO 195 CALL transfert_request(f_du,req_e1_vect) 196 dumax = max_vel(f_du) 197 CALL rescale_field(1./dumax, f_du, f_u) 198 IF (is_master) PRINT *,"gradrot : it :",it ,": dumax",dumax 199 ENDDO 200 343 201 cgradrot=dumax**(-1./nitergrot) 202 IF (is_master) PRINT *,"gradrot : dumax",dumax 344 203 IF (is_master) PRINT *,"cgradrot : ",cgradrot 345 204 346 347 348 !$OMP BARRIER 349 !$OMP MASTER 350 DO ind=1,ndomain 351 CALL swap_dimensions(ind) 352 CALL swap_geometry(ind) 353 theta=f_theta(ind) 354 355 ! set random seed to get reproductibility when using a different number of process 356 CALL RANDOM_SEED(size=M) 357 CALL RANDOM_SEED(put=(/(i,i=1,M)/)) 358 359 DO j=jj_begin,jj_end 360 DO i=ii_begin,ii_end 361 ij=(j-1)*iim+i 362 CALL RANDOM_NUMBER(r) 363 theta(ij)=r-0.5 364 ENDDO 365 ENDDO 366 ENDDO 367 !$OMP END MASTER 368 !$OMP BARRIER 369 370 DO it=1,20 371 372 dthetamax=0 373 DO iter=1,niterdivgrad 374 CALL transfert_request(f_theta,req_i1) 375 DO ind=1,ndomain 376 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 377 CALL swap_dimensions(ind) 378 CALL swap_geometry(ind) 379 theta=f_theta(ind) 380 dtheta=f_dtheta(ind) 381 CALL compute_divgrad(theta,dtheta,1,1) 382 theta=dtheta 383 ENDDO 384 ENDDO 385 386 CALL transfert_request(f_dtheta,req_i1) 387 388 DO ind=1,ndomain 389 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 390 CALL swap_dimensions(ind) 391 CALL swap_geometry(ind) 392 theta=f_theta(ind) 393 dtheta=f_dtheta(ind) 394 395 DO j=jj_begin,jj_end 396 DO i=ii_begin,ii_end 397 ij=(j-1)*iim+i 398 dthetamax=MAX(dthetamax,ABS(dtheta(ij))) 205 !----------------- Compute cdivgrad -------------------- 206 207 cdivgrad=-1. 208 CALL random_scalar(f_theta) 209 DO it=1,20 210 DO iter=1,niterdivgrad 211 CALL transfert_request(f_theta,req_i1) 212 DO ind=1,ndomain 213 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 214 CALL swap_dimensions(ind) 215 CALL swap_geometry(ind) 216 theta=f_theta(ind) 217 dtheta=f_dtheta(ind) 218 CALL compute_divgrad(theta,dtheta,1,1) 399 219 ENDDO 400 ENDDO 401 ENDDO 402 403 IF (using_mpi) THEN 404 CALL reduce_max_omp(dthetamax ,dthetamax1) 405 !$OMP MASTER 406 CALL MPI_ALLREDUCE(dthetamax1,dthetamax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 407 !$OMP END MASTER 408 CALL bcast_omp(dthetamax) 409 ELSE 410 CALL allreduce_max_omp(dthetamax,dthetamax1) 411 dumax=dumax1 412 ENDIF 413 414 IF (is_master) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax 415 416 DO ind=1,ndomain 417 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 418 CALL swap_dimensions(ind) 419 CALL swap_geometry(ind) 420 theta=f_theta(ind) 421 dtheta=f_dtheta(ind) 422 theta=dtheta/dthetamax 423 ENDDO 424 ENDDO 425 220 ENDDO 221 CALL transfert_request(f_dtheta,req_i1) 222 dthetamax = max_scalar(f_dtheta) 223 IF (is_master) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax 224 CALL rescale_field(1./dthetamax, f_dtheta, f_theta) 225 END DO 226 227 cdivgrad=dthetamax**(-1./niterdivgrad) 426 228 IF (is_master) PRINT *,"divgrad : divgrad",dthetamax 427 428 cdivgrad=dthetamax**(-1./niterdivgrad)429 229 IF (is_master) PRINT *,"cdivgrad : ",cdivgrad 430 230 431 231 END SUBROUTINE dissip_constants 232 233 SUBROUTINE dissip_profile 234 USE disvert_mod 235 INTEGER :: l 236 REAL(rstd) :: mintau, zz, zvert, fact 432 237 fact=2 433 238 DO l=1,llm … … 441 246 tau_gradrot(l) = tau_gradrot(l)/zvert 442 247 tau_divgrad(l) = tau_divgrad(l)/zvert 443 END DO444 248 END DO 249 445 250 mintau=tau_graddiv(1) 446 251 DO l=1,llm 447 mintau=MIN(mintau,tau_graddiv(l))448 mintau=MIN(mintau,tau_gradrot(l))449 mintau=MIN(mintau,tau_divgrad(l))450 END DO451 252 mintau=MIN(mintau,tau_graddiv(l)) 253 mintau=MIN(mintau,tau_gradrot(l)) 254 mintau=MIN(mintau,tau_divgrad(l)) 255 END DO 256 452 257 IF(rayleigh_friction_type>0) mintau=MIN(mintau,rayleigh_tau) 453 258 454 259 IF(mintau>0) THEN 455 260 IF (itau_dissip==0) THEN 456 IF (is_master) PRINT *,"init_dissip: Automatic computation of itau_dissip..."457 itau_dissip=INT(mintau/dt)261 IF (is_master) PRINT *,"init_dissip: Automatic computation of itau_dissip..." 262 itau_dissip=INT(mintau/dt) 458 263 ENDIF 459 264 ELSE … … 464 269 dtdissip=itau_dissip*dt 465 270 IF (is_master) THEN 466 PRINT *,"init_dissip: rayleigh_tau",rayleigh_tau, "mintau ",mintau467 PRINT *,"init_dissip: itau_dissip",itau_dissip," dtdissip ",dtdissip271 PRINT *,"init_dissip: rayleigh_tau",rayleigh_tau, "mintau ",mintau 272 PRINT *,"init_dissip: itau_dissip",itau_dissip," dtdissip ",dtdissip 468 273 ENDIF 469 470 END SUBROUTINE init_dissip 471 274 275 END SUBROUTINE dissip_profile 472 276 473 277 SUBROUTINE dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_ue, f_dtheta_rhodz,f_due) 474 USE icosa475 USE theta2theta_rhodz_mod476 USE pression_mod477 USE exner_mod478 USE geopotential_mod479 USE trace480 USE time_mod481 USE omp_para482 IMPLICIT NONE483 278 TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_phis(:), f_geopot(:) 484 279 TYPE(t_field),POINTER :: f_theta_rhodz(:), f_dtheta_rhodz(:) … … 492 287 REAL(rstd),POINTER :: dtheta_rhodz_diss(:,:) 493 288 494 INTEGER :: ind, shear 495 INTEGER :: l,ij,nn 289 INTEGER :: ind, l,ij,nn 496 290 497 291 !$OMP BARRIER … … 500 294 CALL gradiv(f_ue,f_due_diss1) 501 295 CALL gradrot(f_ue,f_due_diss2) 502 503 296 CALL divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz_diss) 504 297 … … 516 309 !DIR$ SIMD 517 310 DO ij=ij_begin,ij_end 518 519 311 due(ij+u_right,l) = -0.5*( due_diss1(ij+u_right,l)/tau_graddiv(l) + due_diss2(ij+u_right,l)/tau_gradrot(l))*itau_dissip 520 312 due(ij+u_lup,l) = -0.5*( due_diss1(ij+u_lup,l) /tau_graddiv(l) + due_diss2(ij+u_lup,l) /tau_gradrot(l))*itau_dissip 521 313 due(ij+u_ldown,l) = -0.5*( due_diss1(ij+u_ldown,l)/tau_graddiv(l) + due_diss2(ij+u_ldown,l)/tau_gradrot(l))*itau_dissip 522 523 314 dtheta_rhodz(ij,l,1) = -0.5*dtheta_rhodz_diss(ij,l)/tau_divgrad(l)*itau_dissip 524 315 ENDDO 525 316 ENDDO 526 527 ! dtheta_rhodz=0 528 ! due=0 529 317 530 318 IF(rayleigh_friction_type>0) THEN 531 IF(rayleigh_friction_type<99) THEN 532 phi=f_geopot(ind) 533 ue=f_ue(ind) 534 DO l=ll_begin,ll_end 319 IF(rayleigh_friction_type<99) THEN 320 phi=f_geopot(ind) 321 ue=f_ue(ind) 322 DO l=ll_begin,ll_end 323 DO ij=ij_begin,ij_end 324 CALL relax(t_right, u_right) 325 CALL relax(t_lup, u_lup) 326 CALL relax(t_ldown, u_ldown) 327 ENDDO 328 END DO 329 ELSE 330 ue=f_ue(ind) 535 331 DO ij=ij_begin,ij_end 536 CALL relax(t_right, u_right) 537 CALL relax(t_lup, u_lup) 538 CALL relax(t_ldown, u_ldown) 332 nn = ij+u_right 333 IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN 334 !print*, "latitude", lat_e(nn)*180./3.14159 335 due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau) 336 ENDIF 337 nn = ij+u_lup 338 IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN 339 due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau) 340 ENDIF 341 nn = ij+u_ldown 342 IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN 343 due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau) 344 ENDIF 539 345 ENDDO 540 END DO 541 ELSE 542 ue=f_ue(ind) 543 DO ij=ij_begin,ij_end 544 nn = ij+u_right 545 IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN 546 !print*, "latitude", lat_e(nn)*180./3.14159 547 due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau) 548 ENDIF 549 nn = ij+u_lup 550 IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN 551 due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau) 552 ENDIF 553 nn = ij+u_ldown 554 IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN 555 due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau) 556 ENDIF 557 ENDDO 558 ENDIF 346 ENDIF 559 347 END IF 560 348 END DO 561 349 562 350 CALL trace_end("dissip") 563 351 564 352 CALL write_dissip_tendencies 565 353 !$OMP BARRIER 566 354 567 CONTAINS 568 569 SUBROUTINE relax(shift_t, shift_u) 570 USE dcmip_initial_conditions_test_1_2_3 571 REAL(rstd) :: z, ulon,ulat, lon,lat, & ! input to test2_schaer_mountain 572 p,hyam,hybm,w,t,phis,ps,rho,q, & ! unused input/output to test2_schaer_mountain 573 fz, u3d(3), uref 574 REAL(rstd), PARAMETER :: zh=2e4,ztop=3e4 ! DCMIP values 575 LOGICAL :: hybrid_eta 576 INTEGER :: shift_u, shift_t, zcoords, nn 577 z = (phi(ij,l)+phi(ij+shift_t,l))/(2.*g) 578 IF(z>zh) THEN ! relax only in the sponge layer z>zh 579 580 nn = ij+shift_u 581 zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm 582 CALL test2_schaer_mountain(lon_e(nn),lat_e(nn),p,z,zcoords,hybrid_eta, & 583 hyam,hybm,shear,ulon,ulat,w,t,phis,ps,rho,q) 584 ! u3d = ulon*elon_e(nn,:) + ulat*elat_e(nn,:) 585 u3d = ulon*elon_e(nn,:) ! ulat=0 586 uref = sum(u3d*ep_e(nn,:)) 587 588 fz = sin((pi/2)*(z-zh)/(ztop-zh)) 589 fz = fz*fz/rayleigh_tau 590 due(nn,l) = due(nn,l) - itau_dissip*fz*(ue(nn,l)-uref) 591 ! fz = 1./rayleigh_tau 592 ! due(nn,l) = due(nn,l) - fz*ue(nn,l) 593 END IF 594 END SUBROUTINE relax 595 596 SUBROUTINE write_dissip_tendencies 597 USE observable_mod, ONLY : f_buf_ulon, f_buf_ulat 598 USE wind_mod 599 USE output_field_mod 600 601 CALL transfert_request(f_due_diss1,req_e1_vect) 602 CALL un2ulonlat(f_due_diss1, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1)))) 603 CALL output_field("dulon_diss1",f_buf_ulon) 604 CALL output_field("dulat_diss1",f_buf_ulat) 605 ! 606 CALL transfert_request(f_due_diss2,req_e1_vect) 607 CALL un2ulonlat(f_due_diss2, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1)))) 608 CALL output_field("dulon_diss2",f_buf_ulon) 609 CALL output_field("dulat_diss2",f_buf_ulat) 610 END SUBROUTINE write_dissip_tendencies 611 612 END SUBROUTINE dissip 613 614 SUBROUTINE gradiv(f_ue,f_due) 615 USE icosa 616 USE trace 617 USE omp_para 618 IMPLICIT NONE 619 TYPE(t_field),POINTER :: f_ue(:) 620 TYPE(t_field),POINTER :: f_due(:) 621 REAL(rstd),POINTER :: ue(:,:) 622 REAL(rstd),POINTER :: due(:,:) 623 INTEGER :: ind 624 INTEGER :: it,l,ij 625 626 CALL trace_start("gradiv") 627 628 DO ind=1,ndomain 355 CONTAINS 356 357 SUBROUTINE relax(shift_t, shift_u) 358 USE dcmip_initial_conditions_test_1_2_3 359 REAL(rstd) :: z, ulon,ulat, lon,lat, & ! input to test2_schaer_mountain 360 p,hyam,hybm,w,t,phis,ps,rho,q, & ! unused input/output to test2_schaer_mountain 361 fz, u3d(3), uref 362 REAL(rstd), PARAMETER :: zh=2e4,ztop=3e4 ! DCMIP values 363 LOGICAL :: hybrid_eta 364 INTEGER :: shift_u, shift_t, zcoords, nn 365 z = (phi(ij,l)+phi(ij+shift_t,l))/(2.*g) 366 IF(z>zh) THEN ! relax only in the sponge layer z>zh 367 nn = ij+shift_u 368 zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm 369 CALL test2_schaer_mountain(lon_e(nn),lat_e(nn),p,z,zcoords,hybrid_eta, & 370 hyam,hybm,rayleigh_shear,ulon,ulat,w,t,phis,ps,rho,q) 371 u3d = ulon*elon_e(nn,:) ! ulat=0 372 uref = sum(u3d*ep_e(nn,:)) 373 374 fz = sin((pi/2)*(z-zh)/(ztop-zh)) 375 fz = fz*fz/rayleigh_tau 376 due(nn,l) = due(nn,l) - itau_dissip*fz*(ue(nn,l)-uref) 377 END IF 378 END SUBROUTINE relax 379 380 SUBROUTINE write_dissip_tendencies 381 USE observable_mod, ONLY : f_buf_ulon, f_buf_ulat 382 USE wind_mod 383 USE output_field_mod 384 385 CALL transfert_request(f_due_diss1,req_e1_vect) 386 CALL un2ulonlat(f_due_diss1, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1)))) 387 CALL output_field("dulon_diss1",f_buf_ulon) 388 CALL output_field("dulat_diss1",f_buf_ulat) 389 ! 390 CALL transfert_request(f_due_diss2,req_e1_vect) 391 CALL un2ulonlat(f_due_diss2, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1)))) 392 CALL output_field("dulon_diss2",f_buf_ulon) 393 CALL output_field("dulat_diss2",f_buf_ulat) 394 END SUBROUTINE write_dissip_tendencies 395 396 END SUBROUTINE dissip 397 398 399 SUBROUTINE gradiv(f_ue,f_due) 400 TYPE(t_field), POINTER :: f_ue(:) 401 TYPE(t_field), POINTER :: f_due(:) 402 REAL(rstd),POINTER :: ue(:,:) 403 REAL(rstd),POINTER :: due(:,:) 404 INTEGER :: ind 405 INTEGER :: it,l,ij 406 407 CALL trace_start("gradiv") 408 409 DO ind=1,ndomain 629 410 IF (.NOT. assigned_domain(ind)) CYCLE 630 411 CALL swap_dimensions(ind) … … 633 414 due=f_due(ind) 634 415 DO l = ll_begin, ll_end 635 !DIR$ SIMD 636 DO ij=ij_begin,ij_end 637 due(ij+u_right,l)=ue(ij+u_right,l) 638 due(ij+u_lup,l)=ue(ij+u_lup,l) 639 due(ij+u_ldown,l)=ue(ij+u_ldown,l) 640 ENDDO 641 ENDDO 642 ENDDO 643 644 DO it=1,nitergdiv 645 416 !DIR$ SIMD 417 DO ij=ij_begin,ij_end 418 due(ij+u_right,l)=ue(ij+u_right,l) 419 due(ij+u_lup,l)=ue(ij+u_lup,l) 420 due(ij+u_ldown,l)=ue(ij+u_ldown,l) 421 ENDDO 422 ENDDO 423 ENDDO 424 425 DO it=1,nitergdiv 646 426 CALL send_message(f_due,req_due) 647 427 CALL wait_message(req_due) 648 649 428 DO ind=1,ndomain 650 IF (.NOT. assigned_domain(ind)) CYCLE651 CALL swap_dimensions(ind)652 CALL swap_geometry(ind)653 due=f_due(ind)654 CALL compute_gradiv(due,due,ll_begin,ll_end)655 ENDDO 656 657 429 IF (.NOT. assigned_domain(ind)) CYCLE 430 CALL swap_dimensions(ind) 431 CALL swap_geometry(ind) 432 due=f_due(ind) 433 CALL compute_gradiv(due,due,ll_begin,ll_end) 434 ENDDO 435 ENDDO 436 658 437 CALL trace_end("gradiv") 659 660 END SUBROUTINE gradiv 661 662 663 SUBROUTINE gradrot(f_ue,f_due) 664 USE icosa 665 USE trace 666 USE omp_para 667 IMPLICIT NONE 668 TYPE(t_field),POINTER :: f_ue(:) 669 TYPE(t_field),POINTER :: f_due(:) 670 REAL(rstd),POINTER :: ue(:,:) 671 REAL(rstd),POINTER :: due(:,:) 672 INTEGER :: ind 673 INTEGER :: it,l,ij 674 438 END SUBROUTINE gradiv 439 440 441 SUBROUTINE gradrot(f_ue,f_due) 442 TYPE(t_field), POINTER :: f_ue(:) 443 TYPE(t_field), POINTER :: f_due(:) 444 REAL(rstd),POINTER :: ue(:,:) 445 REAL(rstd),POINTER :: due(:,:) 446 INTEGER :: ind, it,l,ij 675 447 CALL trace_start("gradrot") 676 448 … … 692 464 693 465 DO it=1,nitergrot 694 695 CALL send_message(f_due,req_due) 696 CALL wait_message(req_due) 697 698 DO ind=1,ndomain 699 IF (.NOT. assigned_domain(ind)) CYCLE 700 CALL swap_dimensions(ind) 701 CALL swap_geometry(ind) 702 due=f_due(ind) 703 CALL compute_gradrot(due,due,ll_begin,ll_end) 704 ENDDO 705 706 ENDDO 707 466 CALL send_message(f_due,req_due) 467 CALL wait_message(req_due) 468 DO ind=1,ndomain 469 IF (.NOT. assigned_domain(ind)) CYCLE 470 CALL swap_dimensions(ind) 471 CALL swap_geometry(ind) 472 due=f_due(ind) 473 CALL compute_gradrot(due,due,ll_begin,ll_end) 474 ENDDO 475 ENDDO 476 708 477 CALL trace_end("gradrot") 709 710 478 END SUBROUTINE gradrot 711 479 712 480 SUBROUTINE divgrad(f_theta,f_dtheta) 713 USE icosa 714 USE trace 715 USE omp_para 716 IMPLICIT NONE 717 TYPE(t_field),POINTER :: f_theta(:) 718 TYPE(t_field),POINTER :: f_dtheta(:) 719 REAL(rstd),POINTER :: theta(:,:) 720 REAL(rstd),POINTER :: dtheta(:,:) 721 INTEGER :: ind 722 INTEGER :: it 481 TYPE(t_field), POINTER :: f_theta(:) 482 TYPE(t_field), POINTER :: f_dtheta(:) 483 REAL(rstd),POINTER :: theta(:,:) 484 REAL(rstd),POINTER :: dtheta(:,:) 485 INTEGER :: ind, it 723 486 724 487 CALL trace_start("divgrad") … … 733 496 ENDDO 734 497 735 DO it=1,niterdivgrad 736 737 CALL transfert_request(f_dtheta,req_i1) 738 498 DO it=1,niterdivgrad 499 CALL transfert_request(f_dtheta,req_i1) 739 500 DO ind=1,ndomain 740 501 IF (.NOT. assigned_domain(ind)) CYCLE … … 744 505 CALL compute_divgrad(dtheta,dtheta,ll_begin,ll_end) 745 506 ENDDO 746 747 507 ENDDO 748 508 749 509 CALL trace_end("divgrad") 750 751 510 END SUBROUTINE divgrad 752 511 753 512 SUBROUTINE divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz) 754 USE icosa 755 USE trace 756 USE omp_para 757 IMPLICIT NONE 758 TYPE(t_field),POINTER :: f_mass(:) 759 TYPE(t_field),POINTER :: f_theta_rhodz(:) 760 TYPE(t_field),POINTER :: f_dtheta_rhodz(:) 513 TYPE(t_field), POINTER :: f_mass(:) 514 TYPE(t_field), POINTER :: f_theta_rhodz(:) 515 TYPE(t_field), POINTER :: f_dtheta_rhodz(:) 761 516 762 517 REAL(rstd),POINTER :: mass(:,:) … … 813 568 ENDDO 814 569 815 816 570 CALL trace_end("divgrad") 817 818 571 END SUBROUTINE divgrad_theta_rhodz 819 572 820 573 SUBROUTINE compute_gradiv(ue,gradivu_e,llb,lle) 821 USE icosa822 IMPLICIT NONE823 574 INTEGER,INTENT(IN) :: llb 824 575 INTEGER,INTENT(IN) :: lle … … 843 594 DO l=llb,lle 844 595 !DIR$ SIMD 845 DO ij=ij_begin,ij_end 846 596 DO ij=ij_begin,ij_end 847 597 gradivu_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*divu_i(ij,l)+ ne(ij+t_right,left)*divu_i(ij+t_right,l) ) 848 849 gradivu_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*divu_i(ij,l)+ ne(ij+t_lup,rdown)*divu_i(ij+t_lup,l)) 850 598 gradivu_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*divu_i(ij,l)+ ne(ij+t_lup,rdown)*divu_i(ij+t_lup,l)) 851 599 gradivu_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*divu_i(ij,l)+ne(ij+t_ldown,rup)*divu_i(ij+t_ldown,l) ) 852 853 600 ENDDO 854 601 ENDDO … … 867 614 868 615 SUBROUTINE compute_divgrad(theta,divgrad_i,llb,lle) 869 USE icosa870 IMPLICIT NONE871 616 INTEGER,INTENT(IN) :: llb 872 617 INTEGER,INTENT(IN) :: lle … … 874 619 REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,llm) 875 620 REAL(rstd) :: grad_e(3*iim*jjm,llb:lle) 876 877 INTEGER :: ij,l 878 879 621 INTEGER :: ij,l 880 622 DO l=llb,lle 881 623 !DIR$ SIMD 882 624 DO ij=ij_begin_ext,ij_end_ext 883 884 grad_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*theta(ij,l)+ ne(ij+t_right,left)*theta(ij+t_right,l) ) 885 886 grad_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*theta(ij,l)+ ne(ij+t_lup,rdown)*theta(ij+t_lup,l )) 887 625 grad_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*theta(ij,l)+ ne(ij+t_right,left)*theta(ij+t_right,l) ) 626 grad_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*theta(ij,l)+ ne(ij+t_lup,rdown)*theta(ij+t_lup,l )) 888 627 grad_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*theta(ij,l)+ne(ij+t_ldown,rup)*theta(ij+t_ldown,l) ) 889 890 ENDDO 891 ENDDO 892 628 ENDDO 629 ENDDO 893 630 894 631 DO l=llb,lle 895 632 !DIR$ SIMD 896 633 DO ij=ij_begin,ij_end 897 898 634 divgrad_i(ij,l)=1./Ai(ij)*(ne(ij,right)*grad_e(ij+u_right,l)*le(ij+u_right) + & 899 635 ne(ij,rup)*grad_e(ij+u_rup,l)*le(ij+u_rup) + & … … 912 648 913 649 END SUBROUTINE compute_divgrad 914 915 650 916 651 SUBROUTINE compute_gradrot(ue,gradrot_e,llb,lle) 917 USE icosa918 IMPLICIT NONE919 652 INTEGER,INTENT(IN) :: llb 920 653 INTEGER,INTENT(IN) :: lle … … 943 676 !DIR$ SIMD 944 677 DO ij=ij_begin,ij_end 945 946 gradrot_e(ij+u_right,l)=1/le(ij+u_right)*ne(ij,right)*(rot_v(ij+z_rdown,l)-rot_v(ij+z_rup,l)) 947 948 gradrot_e(ij+u_lup,l)=1/le(ij+u_lup)*ne(ij,lup)*(rot_v(ij+z_up,l)-rot_v(ij+z_lup,l)) 949 950 gradrot_e(ij+u_ldown,l)=1/le(ij+u_ldown)*ne(ij,ldown)*(rot_v(ij+z_ldown,l)-rot_v(ij+z_down,l)) 951 678 gradrot_e(ij+u_right,l)=1/le(ij+u_right)*ne(ij,right)*(rot_v(ij+z_rdown,l)-rot_v(ij+z_rup,l)) 679 gradrot_e(ij+u_lup,l) =1/le(ij+u_lup) *ne(ij,lup) *(rot_v(ij+z_up,l) -rot_v(ij+z_lup,l)) 680 gradrot_e(ij+u_ldown,l)=1/le(ij+u_ldown)*ne(ij,ldown)*(rot_v(ij+z_ldown,l)-rot_v(ij+z_down,l)) 952 681 ENDDO 953 682 ENDDO … … 955 684 DO l=llb,lle 956 685 !DIR$ SIMD 957 DO ij=ij_begin,ij_end 958 686 DO ij=ij_begin,ij_end 959 687 gradrot_e(ij+u_right,l)=-gradrot_e(ij+u_right,l)*cgradrot 960 688 gradrot_e(ij+u_lup,l)=-gradrot_e(ij+u_lup,l)*cgradrot 961 689 gradrot_e(ij+u_ldown,l)=-gradrot_e(ij+u_ldown,l)*cgradrot 962 963 690 ENDDO 964 691 ENDDO … … 966 693 END SUBROUTINE compute_gradrot 967 694 695 !----------------------- Utility routines ------------------ 696 697 SUBROUTINE global_max(dumax) 698 USE mpi_mod 699 USE mpipara 700 REAL(rstd) :: dumax, dumax1 701 IF (using_mpi) THEN 702 CALL reduce_max_omp(dumax,dumax1) 703 !$OMP MASTER 704 CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 705 !$OMP END MASTER 706 CALL bcast_omp(dumax) 707 ELSE 708 CALL allreduce_max_omp(dumax,dumax1) 709 dumax=dumax1 710 ENDIF 711 END SUBROUTINE global_max 712 713 FUNCTION max_vel(f_du) 714 TYPE(t_field) :: f_du(:) 715 REAL(rstd),POINTER :: du(:) 716 INTEGER :: ind, i,j,ij 717 REAL(rstd) :: max_vel, dumax 718 719 dumax=0. 720 DO ind=1,ndomain 721 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 722 CALL swap_dimensions(ind) 723 CALL swap_geometry(ind) 724 du=f_du(ind) 725 726 DO j=jj_begin,jj_end 727 DO i=ii_begin,ii_end 728 ij=(j-1)*iim+i 729 if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right))) 730 if (le(ij+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_lup))) 731 if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown))) 732 ENDDO 733 ENDDO 734 ENDDO 735 CALL global_max(dumax) 736 max_vel=dumax 737 END FUNCTION max_vel 738 739 FUNCTION max_scalar(f_dtheta) 740 TYPE(t_field) :: f_dtheta(:) 741 REAL(rstd),POINTER :: dtheta(:) 742 INTEGER :: ind, i,j,ij 743 REAL(rstd) :: max_scalar, dthetamax 744 DO ind=1,ndomain 745 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 746 CALL swap_dimensions(ind) 747 dtheta=f_dtheta(ind) 748 DO j=jj_begin,jj_end 749 DO i=ii_begin,ii_end 750 ij=(j-1)*iim+i 751 dthetamax=MAX(dthetamax,ABS(dtheta(ij))) 752 ENDDO 753 ENDDO 754 ENDDO 755 CALL global_max(dthetamax) 756 max_scalar=dthetamax 757 END FUNCTION max_scalar 758 759 SUBROUTINE random_vel(f_u) 760 TYPE(t_field) :: f_u(:) 761 REAL(rstd),POINTER :: u(:) 762 INTEGER :: M, ind, i,j,ij 763 REAL(rstd) :: r 764 765 !$OMP BARRIER 766 !$OMP MASTER 767 DO ind=1,ndomain 768 CALL swap_dimensions(ind) 769 CALL swap_geometry(ind) 770 u=f_u(ind) 771 772 ! set random seed to get reproductibility when using a different number of process 773 CALL RANDOM_SEED(size=M) 774 CALL RANDOM_SEED(put=(/(i,i=1,M)/)) 775 776 DO j=jj_begin,jj_end 777 DO i=ii_begin,ii_end 778 ij=(j-1)*iim+i 779 CALL RANDOM_NUMBER(r) 780 u(ij+u_right)=r-0.5 781 CALL RANDOM_NUMBER(r) 782 u(ij+u_lup)=r-0.5 783 CALL RANDOM_NUMBER(r) 784 u(ij+u_ldown)=r-0.5 785 ENDDO 786 ENDDO 787 ENDDO 788 !$OMP END MASTER 789 !$OMP BARRIER 790 791 END SUBROUTINE random_vel 792 793 SUBROUTINE random_scalar(f_theta) 794 TYPE(t_field) :: f_theta(:) 795 REAL(rstd),POINTER :: theta(:) 796 INTEGER :: M, ind, i,j,ij 797 REAL(rstd) :: r 798 !$OMP BARRIER 799 !$OMP MASTER 800 DO ind=1,ndomain 801 CALL swap_dimensions(ind) 802 CALL swap_geometry(ind) 803 theta=f_theta(ind) 804 ! set random seed to get reproductibility when using a different number of process 805 CALL RANDOM_SEED(size=M) 806 CALL RANDOM_SEED(put=(/(i,i=1,M)/)) 807 808 DO j=jj_begin,jj_end 809 DO i=ii_begin,ii_end 810 ij=(j-1)*iim+i 811 CALL RANDOM_NUMBER(r) 812 theta(ij)=r-0.5 813 ENDDO 814 ENDDO 815 ENDDO 816 !$OMP END MASTER 817 !$OMP BARRIER 818 END SUBROUTINE random_scalar 819 820 SUBROUTINE rescale_field(scal, f_du, f_u) 821 TYPE(t_field) :: f_du(:), f_u(:) 822 REAL(rstd),POINTER :: du(:), u(:) 823 REAL(rstd) :: scal 824 INTEGER :: ind 825 DO ind=1,ndomain 826 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 827 CALL swap_dimensions(ind) 828 u=f_u(ind) 829 du=f_du(ind) 830 u=scal*du 831 ENDDO 832 END SUBROUTINE rescale_field 968 833 969 834 END MODULE dissip_gcm_mod 970
Note: See TracChangeset
for help on using the changeset viewer.