Changeset 360 for codes/icosagcm/trunk/src/timeloop_gcm.f90
- Timestamp:
- 09/08/15 16:23:06 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/timeloop_gcm.f90
r354 r360 1 1 MODULE timeloop_gcm_mod 2 USE transfert_mod3 2 USE icosa 3 USE disvert_mod 4 USE trace 5 USE omp_para 6 USE euler_scheme_mod 7 USE explicit_scheme_mod 8 USE hevi_scheme_mod 4 9 IMPLICIT NONE 5 10 PRIVATE 6 11 12 INTEGER, PARAMETER :: itau_sync=10 13 TYPE(t_message),SAVE :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q0 14 7 15 PUBLIC :: init_timeloop, timeloop 8 9 INTEGER, PARAMETER :: euler=1, rk4=2, mlf=3, rk25=410 INTEGER, PARAMETER :: itau_sync=1011 REAL(rstd), DIMENSION(4), PARAMETER :: coef_rk4 = (/ .25, 1./3., .5, 1. /)12 REAL(rstd), DIMENSION(5), PARAMETER :: coef_rk25 = (/ .25, 1./6., 3./8., .5, 1. /)13 14 TYPE(t_message),SAVE :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q015 16 TYPE(t_field),POINTER,SAVE :: f_geopot(:)17 TYPE(t_field),POINTER,SAVE :: f_q(:)18 TYPE(t_field),POINTER,SAVE :: f_rhodz(:), f_mass(:), f_massm1(:), f_massm2(:), f_dmass(:)19 TYPE(t_field),POINTER,SAVE :: f_phis(:), f_ps(:),f_psm1(:), f_psm2(:), f_dps(:)20 TYPE(t_field),POINTER,SAVE :: f_u(:),f_um1(:),f_um2(:), f_du(:)21 TYPE(t_field),POINTER,SAVE :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:), f_dtheta_rhodz(:)22 TYPE(t_field),POINTER,SAVE :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:)23 24 INTEGER,SAVE :: nb_stage, matsuno_period, scheme25 !$OMP THREADPRIVATE(nb_stage, matsuno_period, scheme)26 16 27 17 CONTAINS 28 18 29 19 SUBROUTINE init_timeloop 30 USE icosa 31 USE dissip_gcm_mod 32 USE observable_mod 33 USE caldyn_mod 34 USE etat0_mod 35 USE disvert_mod 36 USE guided_mod 37 USE advect_tracer_mod 38 USE physics_mod 39 USE mpipara 40 USE omp_para 41 USE trace 42 USE transfert_mod 43 USE check_conserve_mod 44 USE output_field_mod 45 USE write_field_mod 46 USE theta2theta_rhodz_mod 47 USE sponge_mod 48 49 CHARACTER(len=255) :: def 50 51 52 IF (xios_output) itau_out=1 53 IF (.NOT. enable_io) itau_out=HUGE(itau_out) 54 55 ! Time-independant orography 20 USE dissip_gcm_mod 21 USE observable_mod 22 USE caldyn_mod 23 USE etat0_mod 24 USE guided_mod 25 USE advect_tracer_mod 26 USE check_conserve_mod 27 USE output_field_mod 28 USE theta2theta_rhodz_mod 29 USE sponge_mod 30 31 CHARACTER(len=255) :: def 32 33 IF (xios_output) itau_out=1 34 IF (.NOT. enable_io) itau_out=HUGE(itau_out) 35 36 def='RK2.5' 37 CALL getin('time_scheme',def) 38 SELECT CASE (TRIM(def)) 39 CASE('euler') 40 scheme_family=explicit 41 scheme=euler 42 nb_stage=1 43 CASE ('runge_kutta') 44 scheme_family=explicit 45 scheme=rk4 46 nb_stage=4 47 CASE ('RK2.5') 48 scheme_family=explicit 49 scheme=rk25 50 nb_stage=5 51 CASE ('leapfrog_matsuno') 52 scheme_family=explicit 53 scheme=mlf 54 matsuno_period=5 55 CALL getin('matsuno_period',matsuno_period) 56 nb_stage=matsuno_period+1 57 CASE('ARK2.3') 58 scheme_family=hevi 59 scheme=ark23 60 nb_stage=3 61 CALL set_coefs_ark23(dt) 62 CASE('ARK3.3') 63 scheme_family=hevi 64 scheme=ark33 65 nb_stage=3 66 CALL set_coefs_ark33(dt) 67 CASE ('none') 68 nb_stage=0 69 CASE default 70 PRINT*,'Bad selector for variable scheme : <', TRIM(def), & 71 ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>,<RK2.5>,<ARK2.3>' 72 STOP 73 END SELECT 74 75 ! Time-independant orography 56 76 CALL allocate_field(f_phis,field_t,type_real,name='phis') 57 ! Trends 58 CALL allocate_field(f_du,field_u,type_real,llm,name='du') 59 CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm,name='dtheta_rhodz') 60 ! Model state at current time step (RK/MLF/Euler) 77 ! Model state at current time step 78 CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot') 61 79 CALL allocate_field(f_ps,field_t,type_real, name='ps') 62 CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot')63 80 CALL allocate_field(f_mass,field_t,type_real,llm,name='mass') 81 CALL allocate_field(f_rhodz,field_t,type_real,llm,name='rhodz') 82 CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,name='theta_rhodz') 64 83 CALL allocate_field(f_u,field_u,type_real,llm,name='u') 65 CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,name='theta_rhodz')66 ! Model state at previous time step (RK/MLF)67 CALL allocate_field(f_um1,field_u,type_real,llm,name='um1')68 CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm,name='theta_rhodzm1')69 ! Tracers70 84 CALL allocate_field(f_q,field_t,type_real,llm,nqtot,'q') 71 CALL allocate_field(f_rhodz,field_t,type_real,llm,name='rhodz') 72 ! Mass fluxes 85 ! Mass fluxes 73 86 CALL allocate_field(f_hflux,field_u,type_real,llm) ! instantaneous mass fluxes 74 87 CALL allocate_field(f_hfluxt,field_u,type_real,llm) ! mass "fluxes" accumulated in time 75 88 CALL allocate_field(f_wflux,field_t,type_real,llm+1) ! vertical mass fluxes 76 CALL allocate_field(f_dmass,field_t,type_real,llm, name='dmass') 77 78 IF(caldyn_eta == eta_mass) THEN ! eta = mass coordinate (default) 89 CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt') 90 91 SELECT CASE(scheme_family) 92 CASE(explicit) 93 ! Trends 79 94 CALL allocate_field(f_dps,field_t,type_real,name='dps') 95 CALL allocate_field(f_dmass,field_t,type_real,llm, name='dmass') 96 CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm,name='dtheta_rhodz') 97 CALL allocate_field(f_du,field_u,type_real,llm,name='du') 98 ! Model state at previous time step (RK/MLF) 80 99 CALL allocate_field(f_psm1,field_t,type_real,name='psm1') 81 CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt')82 ! the following are unused but must point to something83 ! f_massm1 => f_mass84 ELSE ! eta = Lagrangian vertical coordinate85 100 CALL allocate_field(f_massm1,field_t,type_real,llm, name='massm1') 86 ! the following are unused but must point to something 87 f_wfluxt => f_wflux 88 f_dps => f_phis 89 f_psm1 => f_phis 90 END IF 91 92 def='runge_kutta' 93 CALL getin('scheme',def) 94 95 SELECT CASE (TRIM(def)) 96 CASE('euler') 97 scheme=euler 98 nb_stage=1 99 CASE ('runge_kutta') 100 scheme=rk4 101 nb_stage=4 102 CASE ('RK2.5') 103 scheme=rk25 104 nb_stage=5 105 CASE ('leapfrog_matsuno') 106 scheme=mlf 107 matsuno_period=5 108 CALL getin('matsuno_period',matsuno_period) 109 nb_stage=matsuno_period+1 110 ! Model state 2 time steps ago (MLF) 111 CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm) 112 CALL allocate_field(f_um2,field_u,type_real,llm) 113 IF(caldyn_eta == eta_mass) THEN ! eta = mass coordinate (default) 114 CALL allocate_field(f_psm2,field_t,type_real) 115 ! the following are unused but must point to something 116 f_massm2 => f_mass 117 ELSE ! eta = Lagrangian vertical coordinate 118 CALL allocate_field(f_massm2,field_t,type_real,llm) 119 ! the following are unused but must point to something 120 f_psm2 => f_phis 121 END IF 122 CASE ('none') 123 nb_stage=0 124 125 CASE default 126 PRINT*,'Bad selector for variable scheme : <', TRIM(def), & 127 ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>' 128 STOP 101 CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm,name='theta_rhodzm1') 102 CALL allocate_field(f_um1,field_u,type_real,llm,name='um1') 103 CASE(hevi) 104 ! Trends 105 CALL allocate_fields(nb_stage,f_dps_slow, field_t,type_real,name='dps_slow') 106 CALL allocate_fields(nb_stage,f_dmass_slow, field_t,type_real,llm, name='dmass_slow') 107 CALL allocate_fields(nb_stage,f_dtheta_rhodz_slow, field_t,type_real,llm,name='dtheta_rhodz_fast') 108 CALL allocate_fields(nb_stage,f_du_slow, field_u,type_real,llm,name='du_slow') 109 CALL allocate_fields(nb_stage,f_du_fast, field_u,type_real,llm,name='du_fast') 110 f_dps => f_dps_slow(:,1) 111 f_du => f_du_slow(:,1) 112 f_dtheta_rhodz => f_dtheta_rhodz_slow(:,1) 113 END SELECT 114 115 SELECT CASE(scheme) 116 CASE(mlf) 117 ! Model state 2 time steps ago (MLF) 118 CALL allocate_field(f_psm2,field_t,type_real) 119 CALL allocate_field(f_massm2,field_t,type_real,llm) 120 CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm) 121 CALL allocate_field(f_um2,field_u,type_real,llm) 129 122 END SELECT 130 123 … … 153 146 154 147 SUBROUTINE timeloop 155 USE icosa 156 USE dissip_gcm_mod 157 USE sponge_mod 158 USE disvert_mod 159 USE caldyn_mod 160 USE caldyn_gcm_mod, ONLY : req_ps, req_mass 161 USE observable_mod 162 USE etat0_mod 163 USE guided_mod 164 USE advect_tracer_mod 165 USE physics_mod 166 USE mpipara 167 USE omp_para 168 USE trace 169 USE transfert_mod 170 USE check_conserve_mod 171 USE xios_mod 172 USE output_field_mod 173 USE write_etat0_mod 174 USE checksum_mod 175 IMPLICIT NONE 176 REAL(rstd),POINTER :: q(:,:,:) 177 REAL(rstd),POINTER :: phis(:), ps(:) ,psm1(:), psm2(:), dps(:) 178 REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:), du(:,:) 179 REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), massm1(:,:), massm2(:,:), dmass(:,:) 180 REAL(rstd),POINTER :: theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:), dtheta_rhodz(:,:) 181 REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:) 148 USE dissip_gcm_mod 149 USE sponge_mod 150 USE observable_mod 151 USE etat0_mod 152 USE guided_mod 153 USE caldyn_mod 154 USE advect_tracer_mod 155 USE physics_mod 156 USE mpipara 157 USE transfert_mod 158 USE check_conserve_mod 159 USE xios_mod 160 USE output_field_mod 161 USE write_etat0_mod 162 USE checksum_mod 163 USE explicit_scheme_mod 164 USE hevi_scheme_mod 165 REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), ps(:) 182 166 183 167 INTEGER :: ind 184 INTEGER :: it,i,j, n, stage168 INTEGER :: it,i,j,l,n, stage 185 169 LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time 186 LOGICAL, PARAMETER :: check=.FALSE. 187 INTEGER :: start_clock 188 INTEGER :: stop_clock 189 INTEGER :: rate_clock 190 INTEGER :: l 170 LOGICAL, PARAMETER :: check_rhodz=.FALSE. 171 INTEGER :: start_clock, stop_clock, rate_clock 191 172 LOGICAL,SAVE :: first_physic=.TRUE. 192 !$OMP THREADPRIVATE(first_physic) 193 194 173 !$OMP THREADPRIVATE(first_physic) 174 195 175 CALL switch_omp_distrib_level 196 176 CALL caldyn_BC(f_phis, f_geopot, f_wflux) ! set constant values in first/last interfaces 197 198 !$OMP BARRIER 199 DO ind=1,ndomain 200 IF (.NOT. assigned_domain(ind)) CYCLE 201 CALL swap_dimensions(ind) 202 CALL swap_geometry(ind) 203 rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind) 204 IF(caldyn_eta==eta_mass) THEN 205 CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps 206 ELSE 207 DO l=ll_begin,ll_end 208 rhodz(:,l)=mass(:,l) 209 ENDDO 210 END IF 211 END DO 212 !$OMP BARRIER 213 fluxt_zero=.TRUE. 214 215 !$OMP MASTER 216 CALL SYSTEM_CLOCK(start_clock) 217 !$OMP END MASTER 218 219 CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,itau0) 220 221 CALL trace_on 222 223 DO it=itau0+1,itau0+itaumax 224 225 CALL check_conserve_detailed('detailed_budget 0', & 226 f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 227 228 IF (xios_output) CALL xios_update_calendar(it) 229 IF (it==itau0+1 .OR. MOD(it,itau_sync)==0) THEN 230 CALL send_message(f_ps,req_ps0) 231 CALL wait_message(req_ps0) 232 CALL send_message(f_mass,req_mass0) 233 CALL wait_message(req_mass0) 234 CALL send_message(f_theta_rhodz,req_theta_rhodz0) 235 CALL wait_message(req_theta_rhodz0) 236 CALL send_message(f_u,req_u0) 237 CALL wait_message(req_u0) 238 CALL send_message(f_q,req_q0) 239 CALL wait_message(req_q0) 240 241 ENDIF 242 243 IF (is_master) PRINT *,"It No :",It," t :",dt*It 244 245 IF (mod(it,itau_out)==0 ) THEN 246 CALL update_time_counter(dt*it) 247 CALL write_output_fields_basic(f_ps, f_u, f_q) 248 ! CALL output_field("q",f_q) 249 ENDIF 250 251 CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q) 252 253 DO stage=1,nb_stage 254 ! CALL checksum(f_ps) 255 ! CALL checksum(f_theta_rhodz) 256 ! CALL checksum(f_mass) 257 CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), & 258 f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q, & 259 f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du) 260 ! CALL checksum(f_dps) 261 ! CALL checksum(f_dtheta_rhodz) 262 ! CALL checksum(f_dmass) 263 SELECT CASE (scheme) 264 CASE(euler) 265 CALL euler_scheme(.TRUE.) 266 CASE (rk4) 267 CALL rk_scheme(stage, coef_rk4) 268 CASE (rk25) 269 CALL rk_scheme(stage, coef_rk25) 270 CASE (mlf) 271 CALL leapfrog_matsuno_scheme(stage) 272 CASE DEFAULT 273 STOP 274 END SELECT 275 END DO 276 277 CALL check_conserve_detailed('detailed_budget 1', & 278 f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 279 280 IF (MOD(it,itau_dissip)==0) THEN 281 282 IF(caldyn_eta==eta_mass) THEN 283 !ym flush ps 284 !$OMP BARRIER 285 DO ind=1,ndomain 286 IF (.NOT. assigned_domain(ind)) CYCLE 287 CALL swap_dimensions(ind) 288 CALL swap_geometry(ind) 289 mass=f_mass(ind); ps=f_ps(ind); 290 CALL compute_rhodz(.TRUE., ps, mass) 291 END DO 292 ENDIF 293 294 CALL dissip(f_u,f_du,f_mass,f_phis, f_theta_rhodz,f_dtheta_rhodz) 295 296 CALL euler_scheme(.FALSE.) ! update only u, theta 297 IF (iflag_sponge > 0) THEN 298 CALL sponge(f_u,f_du,f_theta_rhodz,f_dtheta_rhodz) 299 CALL euler_scheme(.FALSE.) ! update only u, theta 300 ENDIF 301 END IF 302 303 CALL check_conserve_detailed('detailed_budget 2', & 304 f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 305 306 IF(MOD(it,itau_adv)==0) THEN 307 308 CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz) ! update q and rhodz after RK step 309 fluxt_zero=.TRUE. 310 311 ! FIXME : check that rhodz is consistent with ps 312 IF (check) THEN 313 DO ind=1,ndomain 314 IF (.NOT. assigned_domain(ind)) CYCLE 315 CALL swap_dimensions(ind) 316 CALL swap_geometry(ind) 317 rhodz=f_rhodz(ind); ps=f_ps(ind); 318 CALL compute_rhodz(.FALSE., ps, rhodz) 319 END DO 320 ENDIF 321 322 END IF 323 324 325 IF (MOD(it,itau_check_conserv)==0) THEN 326 CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 327 ENDIF 328 329 IF (MOD(it,itau_physics)==0) THEN 330 CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q) 331 332 !$OMP MASTER 333 IF (first_physic) CALL SYSTEM_CLOCK(start_clock) 334 !$OMP END MASTER 335 first_physic=.FALSE. 336 ENDIF 337 338 ENDDO 339 340 CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q) 341 342 CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 343 344 !$OMP MASTER 345 CALL SYSTEM_CLOCK(stop_clock) 346 CALL SYSTEM_CLOCK(count_rate=rate_clock) 347 348 IF (mpi_rank==0) THEN 349 PRINT *,"Time elapsed : ",(stop_clock-start_clock)*1./rate_clock 350 ENDIF 351 !$OMP END MASTER 352 353 CONTAINS 354 355 SUBROUTINE Euler_scheme(with_dps) 356 IMPLICIT NONE 357 LOGICAL :: with_dps 358 INTEGER :: ind 359 INTEGER :: i,j,ij,l 360 CALL trace_start("Euler_scheme") 361 177 178 !$OMP BARRIER 362 179 DO ind=1,ndomain 363 180 IF (.NOT. assigned_domain(ind)) CYCLE 364 181 CALL swap_dimensions(ind) 365 182 CALL swap_geometry(ind) 366 367 IF(with_dps) THEN ! update ps/mass 368 IF(caldyn_eta==eta_mass) THEN ! update ps 369 ps=f_ps(ind) ; dps=f_dps(ind) ; 370 IF (is_omp_first_level) THEN 371 !$SIMD 372 DO ij=ij_begin,ij_end 373 ps(ij)=ps(ij)+dt*dps(ij) 374 ENDDO 375 ENDIF 376 ELSE ! update mass 377 mass=f_mass(ind) ; dmass=f_dmass(ind) ; 378 DO l=ll_begin,ll_end 379 !$SIMD 380 DO ij=ij_begin,ij_end 381 mass(ij,l)=mass(ij,l)+dt*dmass(ij,l) 382 ENDDO 383 END DO 384 END IF 385 386 hflux=f_hflux(ind); hfluxt=f_hfluxt(ind) 387 wflux=f_wflux(ind); wfluxt=f_wfluxt(ind) 388 CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind)) 389 END IF ! update ps/mass 390 391 u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) 392 du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind) 393 394 DO l=ll_begin,ll_end 395 !$SIMD 396 DO ij=ij_begin,ij_end 397 u(ij+u_right,l)=u(ij+u_right,l)+dt*du(ij+u_right,l) 398 u(ij+u_lup,l)=u(ij+u_lup,l)+dt*du(ij+u_lup,l) 399 u(ij+u_ldown,l)=u(ij+u_ldown,l)+dt*du(ij+u_ldown,l) 400 theta_rhodz(ij,l)=theta_rhodz(ij,l)+dt*dtheta_rhodz(ij,l) 401 ENDDO 402 ENDDO 403 ENDDO 404 405 CALL trace_end("Euler_scheme") 406 407 END SUBROUTINE Euler_scheme 408 409 SUBROUTINE RK_scheme(stage,coef) 410 IMPLICIT NONE 411 INTEGER :: ind, stage 412 REAL(rstd), INTENT(IN) :: coef(:) 413 REAL(rstd) :: tau 414 INTEGER :: i,j,ij,l 415 416 CALL trace_start("RK_scheme") 417 418 tau = dt*coef(stage) 419 420 ! if mass coordinate, deal with ps first on one core 421 IF(caldyn_eta==eta_mass) THEN 422 IF (is_omp_first_level) THEN 423 424 DO ind=1,ndomain 425 IF (.NOT. assigned_domain(ind)) CYCLE 426 CALL swap_dimensions(ind) 427 CALL swap_geometry(ind) 428 ps=f_ps(ind) ; psm1=f_psm1(ind) ; dps=f_dps(ind) 429 430 IF (stage==1) THEN ! first stage : save model state in XXm1 431 !$SIMD 432 DO ij=ij_begin,ij_end 433 psm1(ij)=ps(ij) 434 ENDDO 435 ENDIF 436 437 ! updates are of the form x1 := x0 + tau*f(x1) 438 !$SIMD 439 DO ij=ij_begin,ij_end 440 ps(ij)=psm1(ij)+tau*dps(ij) 441 ENDDO 442 ENDDO 443 ENDIF 444 ! CALL send_message(f_ps,req_ps) 445 !ym no overlap for now 446 ! CALL wait_message(req_ps) 447 448 ELSE ! Lagrangian coordinate, deal with mass 449 DO ind=1,ndomain 450 IF (.NOT. assigned_domain(ind)) CYCLE 451 CALL swap_dimensions(ind) 452 CALL swap_geometry(ind) 453 mass=f_mass(ind); dmass=f_dmass(ind); massm1=f_massm1(ind) 454 455 IF (stage==1) THEN ! first stage : save model state in XXm1 456 DO l=ll_begin,ll_end 457 !$SIMD 458 DO ij=ij_begin,ij_end 459 massm1(ij,l)=mass(ij,l) 460 ENDDO 461 ENDDO 462 END IF 463 464 ! updates are of the form x1 := x0 + tau*f(x1) 465 DO l=ll_begin,ll_end 466 !$SIMD 467 DO ij=ij_begin,ij_end 468 mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l) 469 ENDDO 470 ENDDO 471 END DO 472 ! CALL send_message(f_mass,req_mass) 473 !ym no overlap for now 474 ! CALL wait_message(req_mass) 475 476 END IF 477 478 ! now deal with other prognostic variables 479 DO ind=1,ndomain 480 IF (.NOT. assigned_domain(ind)) CYCLE 481 CALL swap_dimensions(ind) 482 CALL swap_geometry(ind) 483 u=f_u(ind) ; du=f_du(ind) ; um1=f_um1(ind) 484 theta_rhodz=f_theta_rhodz(ind) 485 theta_rhodzm1=f_theta_rhodzm1(ind) 486 dtheta_rhodz=f_dtheta_rhodz(ind) 487 488 IF (stage==1) THEN ! first stage : save model state in XXm1 489 DO l=ll_begin,ll_end 490 !$SIMD 491 DO ij=ij_begin,ij_end 492 um1(ij+u_right,l)=u(ij+u_right,l) 493 um1(ij+u_lup,l)=u(ij+u_lup,l) 494 um1(ij+u_ldown,l)=u(ij+u_ldown,l) 495 theta_rhodzm1(ij,l)=theta_rhodz(ij,l) 496 ENDDO 497 ENDDO 498 END IF 499 500 DO l=ll_begin,ll_end 501 !$SIMD 502 DO ij=ij_begin,ij_end 503 u(ij+u_right,l)=um1(ij+u_right,l)+tau*du(ij+u_right,l) 504 u(ij+u_lup,l)=um1(ij+u_lup,l)+tau*du(ij+u_lup,l) 505 u(ij+u_ldown,l)=um1(ij+u_ldown,l)+tau*du(ij+u_ldown,l) 506 theta_rhodz(ij,l)=theta_rhodzm1(ij,l)+tau*dtheta_rhodz(ij,l) 507 ENDDO 508 ENDDO 509 510 IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage 511 hflux=f_hflux(ind); hfluxt=f_hfluxt(ind) 512 wflux=f_wflux(ind); wfluxt=f_wfluxt(ind) 513 CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind)) 514 END IF 515 END DO 516 517 CALL trace_end("RK_scheme") 518 519 END SUBROUTINE RK_scheme 520 521 SUBROUTINE leapfrog_matsuno_scheme(stage) 522 IMPLICIT NONE 523 INTEGER :: ind, stage 524 REAL :: tau 525 526 CALL trace_start("leapfrog_matsuno_scheme") 527 528 tau = dt/nb_stage 529 DO ind=1,ndomain 530 IF (.NOT. assigned_domain(ind)) CYCLE 531 CALL swap_dimensions(ind) 532 CALL swap_geometry(ind) 533 534 ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) 535 psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind) 536 psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind) 537 dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind) 538 539 540 IF (stage==1) THEN 541 psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) 542 ps(:)=psm1(:)+tau*dps(:) 543 u(:,:)=um1(:,:)+tau*du(:,:) 544 theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:) 545 546 ELSE IF (stage==2) THEN 547 548 ps(:)=psm1(:)+tau*dps(:) 549 u(:,:)=um1(:,:)+tau*du(:,:) 550 theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:) 551 552 psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 553 psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 554 555 ELSE 556 557 ps(:)=psm2(:)+2*tau*dps(:) 558 u(:,:)=um2(:,:)+2*tau*du(:,:) 559 theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*tau*dtheta_rhodz(:,:) 560 561 psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 562 psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 563 564 ENDIF 565 566 ENDDO 567 CALL trace_end("leapfrog_matsuno_scheme") 568 569 END SUBROUTINE leapfrog_matsuno_scheme 570 571 END SUBROUTINE timeloop 572 573 SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero) 574 USE icosa 575 USE omp_para 576 USE disvert_mod 577 IMPLICIT NONE 578 REAL(rstd), INTENT(IN) :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1) 579 REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1) 580 REAL(rstd), INTENT(IN) :: tau 581 LOGICAL, INTENT(INOUT) :: fluxt_zero 582 INTEGER :: l,i,j,ij 583 584 IF(fluxt_zero) THEN 585 586 fluxt_zero=.FALSE. 587 588 DO l=ll_begin,ll_end 589 !$SIMD 590 DO ij=ij_begin_ext,ij_end_ext 591 hfluxt(ij+u_right,l) = tau*hflux(ij+u_right,l) 592 hfluxt(ij+u_lup,l) = tau*hflux(ij+u_lup,l) 593 hfluxt(ij+u_ldown,l) = tau*hflux(ij+u_ldown,l) 594 ENDDO 595 ENDDO 596 597 IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag 598 DO l=ll_begin,ll_endp1 599 !$SIMD 600 DO ij=ij_begin,ij_end 601 wfluxt(ij,l) = tau*wflux(ij,l) 602 ENDDO 183 rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind) 184 IF(caldyn_eta==eta_mass) THEN 185 CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps 186 ELSE 187 DO l=ll_begin,ll_end 188 rhodz(:,l)=mass(:,l) 603 189 ENDDO 604 190 END IF 605 606 ELSE 607 608 DO l=ll_begin,ll_end 609 !$SIMD 610 DO ij=ij_begin_ext,ij_end_ext 611 hfluxt(ij+u_right,l) = hfluxt(ij+u_right,l)+tau*hflux(ij+u_right,l) 612 hfluxt(ij+u_lup,l) = hfluxt(ij+u_lup,l)+tau*hflux(ij+u_lup,l) 613 hfluxt(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)+tau*hflux(ij+u_ldown,l) 614 ENDDO 615 ENDDO 616 617 IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag 618 DO l=ll_begin,ll_endp1 619 !$SIMD 620 DO ij=ij_begin,ij_end 621 wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l) 622 ENDDO 623 ENDDO 191 END DO 192 !$OMP BARRIER 193 fluxt_zero=.TRUE. 194 195 !$OMP MASTER 196 CALL SYSTEM_CLOCK(start_clock) 197 !$OMP END MASTER 198 199 CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,itau0) 200 201 CALL trace_on 202 203 DO it=itau0+1,itau0+itaumax 204 205 CALL check_conserve_detailed('detailed_budget 0', & 206 f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 207 208 IF (xios_output) CALL xios_update_calendar(it) 209 IF (it==itau0+1 .OR. MOD(it,itau_sync)==0) THEN 210 CALL send_message(f_ps,req_ps0) 211 CALL wait_message(req_ps0) 212 CALL send_message(f_mass,req_mass0) 213 CALL wait_message(req_mass0) 214 CALL send_message(f_theta_rhodz,req_theta_rhodz0) 215 CALL wait_message(req_theta_rhodz0) 216 CALL send_message(f_u,req_u0) 217 CALL wait_message(req_u0) 218 CALL send_message(f_q,req_q0) 219 CALL wait_message(req_q0) 220 ENDIF 221 222 IF (is_master) PRINT *,"It No :",It," t :",dt*It 223 224 IF (mod(it,itau_out)==0 ) THEN 225 CALL update_time_counter(dt*it) 226 CALL write_output_fields_basic(f_ps, f_u, f_q) 227 ENDIF 228 229 CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q) 230 231 SELECT CASE(scheme_family) 232 CASE(explicit) 233 CALL explicit_scheme(it, fluxt_zero) 234 CASE(hevi) 235 CALL HEVI_scheme(it, fluxt_zero) 236 END SELECT 237 238 CALL check_conserve_detailed('detailed_budget 1', & 239 f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 240 241 IF (MOD(it,itau_dissip)==0) THEN 242 243 IF(caldyn_eta==eta_mass) THEN 244 !ym flush ps 245 !$OMP BARRIER 246 DO ind=1,ndomain 247 IF (.NOT. assigned_domain(ind)) CYCLE 248 CALL swap_dimensions(ind) 249 CALL swap_geometry(ind) 250 mass=f_mass(ind); ps=f_ps(ind); 251 CALL compute_rhodz(.TRUE., ps, mass) 252 END DO 253 ENDIF 254 255 CALL dissip(f_u,f_du,f_mass,f_phis, f_theta_rhodz,f_dtheta_rhodz) 256 257 CALL euler_scheme(.FALSE.) ! update only u, theta 258 IF (iflag_sponge > 0) THEN 259 CALL sponge(f_u,f_du,f_theta_rhodz,f_dtheta_rhodz) 260 CALL euler_scheme(.FALSE.) ! update only u, theta 261 END IF 624 262 END IF 625 626 END IF 627 628 END SUBROUTINE accumulate_fluxes 263 264 CALL check_conserve_detailed('detailed_budget 2', & 265 f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 266 267 IF(MOD(it,itau_adv)==0) THEN 268 CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz) ! update q and rhodz after RK step 269 fluxt_zero=.TRUE. 270 ! FIXME : check that rhodz is consistent with ps 271 IF (check_rhodz) THEN 272 DO ind=1,ndomain 273 IF (.NOT. assigned_domain(ind)) CYCLE 274 CALL swap_dimensions(ind) 275 CALL swap_geometry(ind) 276 rhodz=f_rhodz(ind); ps=f_ps(ind); 277 CALL compute_rhodz(.FALSE., ps, rhodz) 278 END DO 279 ENDIF 280 281 END IF 282 283 IF (MOD(it,itau_check_conserv)==0) THEN 284 CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 285 ENDIF 286 287 IF (MOD(it,itau_physics)==0) THEN 288 CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q) 289 !$OMP MASTER 290 IF (first_physic) CALL SYSTEM_CLOCK(start_clock) 291 !$OMP END MASTER 292 first_physic=.FALSE. 293 END IF 294 295 END DO 296 297 CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q) 298 299 CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 300 301 !$OMP MASTER 302 CALL SYSTEM_CLOCK(stop_clock) 303 CALL SYSTEM_CLOCK(count_rate=rate_clock) 304 305 IF (mpi_rank==0) THEN 306 PRINT *,"Time elapsed : ",(stop_clock-start_clock)*1./rate_clock 307 ENDIF 308 !$OMP END MASTER 309 310 ! CONTAINS 311 END SUBROUTINE timeloop 629 312 630 313 END MODULE timeloop_gcm_mod
Note: See TracChangeset
for help on using the changeset viewer.