Changeset 714 for codes/icosagcm
- Timestamp:
- 08/03/18 16:53:37 (6 years ago)
- Location:
- codes/icosagcm/devel/src
- Files:
-
- 1 added
- 11 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/base/earth_const.f90
r628 r714 13 13 REAL(rstd),SAVE :: Treff=273. 14 14 REAL(rstd),SAVE :: preff=101325. 15 REAL(rstd),SAVE :: pa=50000. 15 REAL(rstd),SAVE :: pa=50000. ! default value set to preff/2 by disvert_std 16 16 REAL(rstd),SAVE :: scale_height=8000. ! atmospheric scale height (m) 17 17 REAL(rstd),SAVE :: scale_factor=1. -
codes/icosagcm/devel/src/diagnostics/check_conserve.f90
r609 r714 14 14 15 15 REAL(rstd),SAVE :: AAM_mass, AAM_mass_old, AAM_vel, AAM_vel_old, & 16 AAM_velp, AAM_velp_old, AAM_velm, AAM_velm_old, & 16 17 AAM_mass_source(3), AAM_vel_source(3) ! read/written only IF is_master 18 REAL(rstd),SAVE :: AAM_vel_plus_source(3), AAM_vel_minus_source(3) 17 19 REAL(rstd),SAVE :: mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0 18 20 !$OMP THREADPRIVATE(check_type, mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0) … … 66 68 INTEGER :: ind,ierr 67 69 REAL(rstd) :: mtot, angtot, rmsdpdt 68 REAL(rstd) :: etot, stot, ang_mass, ang_vel, rmsvtot, ztot70 REAL(rstd) :: etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsvtot, ztot 69 71 70 72 CALL transfert_request(f_ue,req_e1_vect) … … 84 86 CALL check_PV(ztot) 85 87 CALL exner(f_ps,f_p,f_pks,f_pk) 86 CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, rmsvtot)88 CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsvtot) 87 89 88 90 IF (is_master) THEN ! is_master = is_omp_master && is_mpi_master … … 90 92 AAM_mass = ang_mass 91 93 AAM_vel = ang_vel 94 AAM_velp = ang_velp 95 AAM_velm = ang_velm 92 96 angtot = ang_mass + ang_vel 93 97 IF ( it == itau0 ) THEN … … 99 103 AAM_mass_old=AAM_mass 100 104 AAM_vel_old=AAM_vel 105 AAM_velp_old=AAM_velp 106 AAM_velm_old=AAM_velm 101 107 AAM_mass_source=0. 102 108 AAM_vel_source=0. 109 AAM_vel_plus_source=0. 110 AAM_vel_minus_source=0. 103 111 END IF 104 112 rmsvtot=SQRT(rmsvtot/mtot) … … 127 135 WRITE(*,'(A,6E12.4)') 'AAM_vel : time,old,new,dissip,dyn,phys', dt*it, AAM_vel_old, AAM_vel, & 128 136 AAM_vel_source(AAM_dissip), AAM_vel_source(AAM_dyn), AAM_vel_source(AAM_phys) 129 WRITE(*,'(A,6E12.4)') 'AAM_ tot : time,old,new,dissip,dyn,phys', dt*it, &130 AAM_vel_ old+AAM_mass_old, AAM_vel+AAM_mass, &131 AAM_vel_source(AAM_dissip)+AAM_mass_source(AAM_dissip), &132 AAM_vel_ source(AAM_dyn) +AAM_mass_source(AAM_dyn), &133 AAM_vel_source(AAM_phys) +AAM_mass_source(AAM_phys)137 WRITE(*,'(A,6E12.4)') 'AAM_vel+ : time,old,new,dissip,dyn,phys', dt*it, AAM_velp_old, AAM_velp, & 138 AAM_vel_plus_source(AAM_dissip), AAM_vel_plus_source(AAM_dyn), AAM_vel_plus_source(AAM_phys) 139 WRITE(*,'(A,6E12.4)') 'AAM_vel- : time,old,new,dissip,dyn,phys', dt*it, AAM_velm_old, AAM_velm, & 140 AAM_vel_minus_source(AAM_dissip), AAM_vel_minus_source(AAM_dyn), AAM_vel_minus_source(AAM_phys) 141 WRITE(*,'(A,6E12.4)') 'AAM_dyn budget mass + vel:', AAM_mass_source(AAM_dyn) + AAM_vel_source(AAM_dyn) 134 142 AAM_mass_old=AAM_mass 135 143 AAM_vel_old=AAM_vel 144 AAM_velp_old=AAM_velp 145 AAM_velm_old=AAM_velm 136 146 AAM_mass_source=0. 137 147 AAM_vel_source=0. 148 AAM_vel_minus_source=0. 149 AAM_vel_plus_source=0. 138 150 END IF 139 151 END IF … … 160 172 REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 161 173 INTEGER::ind,ierr 162 REAL(rstd) :: mtot, ztot, rmsdpdt, etot,stot,rmsv, ang_mass, ang_vel 174 REAL(rstd) :: mtot, ztot, rmsdpdt, etot,stot,rmsv, ang_mass, ang_vel, ang_velp, ang_velm 163 175 164 176 IF(check_type == check_detailed) THEN … … 175 187 END DO 176 188 CALL exner(f_ps,f_p,f_pks,f_pk) 177 CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, rmsv)189 CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsv) 178 190 179 191 IF (is_master) THEN … … 186 198 AAM_mass_source(tag) = AAM_mass_source(tag) + ang_mass - AAM_mass 187 199 AAM_vel_source(tag) = AAM_vel_source(tag) + ang_vel - AAM_vel 200 AAM_vel_plus_source(tag) = AAM_vel_plus_source(tag) + ang_velp - AAM_velp 201 AAM_vel_minus_source(tag) = AAM_vel_minus_source(tag) + ang_velm - AAM_velm 202 !if ((ang_vel - AAM_vel) .gt. 0.) then 203 ! AAM_vel_plus_source(tag) = AAM_vel_plus_source(tag) + ang_vel - AAM_vel 204 !else 205 ! AAM_vel_minus_source(tag) = AAM_vel_minus_source(tag) + ang_vel - AAM_vel 206 !endif 188 207 AAM_mass = ang_mass 189 208 AAM_vel = ang_vel 209 AAM_velp = ang_velp 210 AAM_velm = ang_velm 190 211 END IF 191 212 END IF … … 252 273 253 274 SUBROUTINE check_energy(f_ue,f_theta_rhodz,f_phis, etot, & 254 stot, AAM_mass_tot, AAM_vel_tot, rmsvtot)275 stot, AAM_mass_tot, AAM_vel_tot, AAM_velp_tot, AAM_velm_tot, rmsvtot) 255 276 USE pression_mod 256 277 USE vorticity_mod … … 258 279 TYPE(t_field), POINTER :: f_theta_rhodz(:) 259 280 TYPE(t_field), POINTER :: f_phis(:) 260 REAL(rstd),INTENT(OUT) :: etot, stot, AAM_mass_tot, AAM_vel_tot, rmsvtot281 REAL(rstd),INTENT(OUT) :: etot, stot, AAM_mass_tot, AAM_vel_tot, AAM_velp_tot, AAM_velm_tot, rmsvtot 261 282 REAL(rstd), POINTER :: ue(:,:) 262 283 REAL(rstd), POINTER :: p(:,:) … … 265 286 REAL(rstd), POINTER :: phis(:) 266 287 REAL(rstd), POINTER :: rhodz(:,:) 267 REAL(rstd) :: e, s, AAM_mass, AAM_vel, rmsv288 REAL(rstd) :: e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv 268 289 INTEGER :: ind 269 290 270 291 e = 0 ; s = 0 ; AAM_mass = 0 ; AAM_vel=0 ; rmsv = 0 292 AAM_velp = 0 ; AAM_velm = 0 271 293 272 294 DO ind=1,ndomain … … 281 303 phis=f_phis(ind) 282 304 CALL compute_energy(ind,ue,p,rhodz,theta_rhodz(:,:,1),pk,phis, & 283 e, s, AAM_mass, AAM_vel, rmsv)305 e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv) 284 306 END DO 285 307 … … 289 311 CALL global_sum(AAM_mass, AAM_mass_tot) 290 312 CALL global_sum(AAM_vel, AAM_vel_tot) 313 CALL global_sum(AAM_velp, AAM_velp_tot) 314 CALL global_sum(AAM_velm, AAM_velm_tot) 291 315 END SUBROUTINE check_energy 292 316 293 SUBROUTINE compute_energy(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, AAM_mass, AAM_vel, rmsv)317 SUBROUTINE compute_energy(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv) 294 318 USE disvert_mod 295 319 USE wind_mod … … 302 326 REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm) 303 327 REAL(rstd),INTENT(IN) :: phis(iim*jjm) 304 REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, rmsv328 REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv 305 329 306 330 REAL(rstd) :: ucenter(iim*jjm,llm,3) … … 335 359 AAM_mass = AAM_mass + rad*cos(lat)*mass_ij*radomeg*cos(lat) 336 360 AAM_vel = AAM_vel + rad*cos(lat)*mass_ij*ulon(ij,l) 361 if (ulon(ij,l) > 0.) then 362 AAM_velp = AAM_velp + rad*cos(lat)*mass_ij*ulon(ij,l) 363 else 364 AAM_velm = AAM_velm + rad*cos(lat)*mass_ij*ulon(ij,l) 365 endif 337 366 END IF 338 367 END DO -
codes/icosagcm/devel/src/diagnostics/observable.f90
r603 r714 15 15 TYPE(t_field),POINTER, SAVE :: f_theta(:) 16 16 17 PUBLIC init_observable, write_output_fields_basic, f_theta, f_buf_i 17 PUBLIC init_observable, write_output_fields_basic, & 18 f_theta, f_buf_i, f_buf_ulon, f_buf_ulat 18 19 19 20 CONTAINS -
codes/icosagcm/devel/src/diagnostics/wind.F90
r612 r714 12 12 CONTAINS 13 13 14 SUBROUTINE un2ulonlat(f_u, f_ulon, f_ulat )14 SUBROUTINE un2ulonlat(f_u, f_ulon, f_ulat, scale_) 15 15 TYPE(t_field), POINTER :: f_u(:) ! IN : normal velocity components on edges 16 16 TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons 17 REAL(rstd),POINTER :: u(:,:), ulon(:,:), ulat(:,:) 17 REAL(rstd),POINTER :: u(:,:), ulon(:,:), ulat(:,:) 18 REAL(rstd), OPTIONAL :: scale_ 19 REAL(rstd) :: scale 18 20 INTEGER :: ind 19 21 scale = MERGE(scale_, 1., PRESENT(scale_)) 20 22 DO ind=1,ndomain 21 23 IF (.NOT. assigned_domain(ind)) CYCLE … … 25 27 ulon=f_ulon(ind) 26 28 ulat=f_ulat(ind) 27 CALL compute_un2ulonlat(u,ulon, ulat )29 CALL compute_un2ulonlat(u,ulon, ulat, scale) 28 30 END DO 29 31 … … 49 51 END SUBROUTINE ulonlat2un 50 52 51 SUBROUTINE compute_wind_centered(ue,ucenter )53 SUBROUTINE compute_wind_centered(ue,ucenter,scale_) 52 54 REAL(rstd) :: ue(3*iim*jjm,llm) 53 55 REAL(rstd) :: ucenter(iim*jjm,llm,3) 56 REAL(rstd), INTENT(IN), OPTIONAL :: scale_ 54 57 INTEGER :: ij,l 55 REAL(rstd) , PARAMETER :: scale=1.56 REAL(rstd) :: fac, ue_le, cx,cy,cz, ux,uy,uz58 REAL(rstd) :: scale,fac, ue_le, cx,cy,cz, ux,uy,uz 59 scale = MERGE(scale_, 1., PRESENT(scale_)) 57 60 #include "../kernels_hex/wind_centered.k90" 58 61 END SUBROUTINE compute_wind_centered … … 328 331 329 332 330 SUBROUTINE compute_un2ulonlat(un, ulon, ulat )333 SUBROUTINE compute_un2ulonlat(un, ulon, ulat, scale) 331 334 REAL(rstd),INTENT(IN) :: un(3*iim*jjm,llm) 332 335 REAL(rstd),INTENT(OUT) :: ulon(iim*jjm,llm) 333 336 REAL(rstd),INTENT(OUT) :: ulat(iim*jjm,llm) 334 335 337 REAL(rstd) :: uc(iim*jjm,llm,3) 336 337 CALL compute_wind_centered(un,uc )338 REAL(rstd) :: scale 339 CALL compute_wind_centered(un,uc,scale) 338 340 CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat) 339 341 -
codes/icosagcm/devel/src/dissip/dissip_gcm.f90
r550 r714 35 35 !$OMP THREADPRIVATE(rayleigh_friction_type) 36 36 !$OMP THREADPRIVATE(rayleigh_shear) 37 REAL, SAVE :: rayleigh_tau 37 REAL, SAVE :: rayleigh_tau, rayleigh_limlat 38 38 !$OMP THREADPRIVATE(rayleigh_tau) 39 !$OMP THREADPRIVATE(rayleigh_limlat) 39 40 40 41 REAL,SAVE :: dtdissip … … 103 104 rayleigh_friction_type=2 104 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' 105 109 CASE DEFAULT 106 110 IF (is_master) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), & … … 117 121 STOP 118 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 119 128 END IF 120 129 … … 444 453 445 454 IF(mintau>0) THEN 446 itau_dissip=INT(mintau/dt) 447 dtdissip=itau_dissip*dt 455 IF (itau_dissip==0) THEN 456 IF (is_master) PRINT *,"init_dissip: Automatic computation of itau_dissip..." 457 itau_dissip=INT(mintau/dt) 458 ENDIF 448 459 ELSE 449 IF (is_master) PRINT *," No dissipation time set, setting itau_dissip to 1000000000"460 IF (is_master) PRINT *,"init_dissip: No dissipation time set, setting itau_dissip to 1000000000" 450 461 itau_dissip=100000000 451 462 END IF 452 463 itau_dissip=MAX(1,itau_dissip) 453 IF (is_master) PRINT *,"rayleigh_tau",rayleigh_tau, "mintau ",mintau, & 454 "itau_dissip",itau_dissip," dtdissip ",dtdissip 464 dtdissip=itau_dissip*dt 465 IF (is_master) THEN 466 PRINT *,"init_dissip: rayleigh_tau",rayleigh_tau, "mintau ",mintau 467 PRINT *,"init_dissip: itau_dissip",itau_dissip," dtdissip ",dtdissip 468 ENDIF 455 469 456 470 END SUBROUTINE init_dissip … … 479 493 480 494 INTEGER :: ind, shear 481 INTEGER :: l,ij 495 INTEGER :: l,ij,nn 482 496 483 497 !$OMP BARRIER … … 515 529 516 530 IF(rayleigh_friction_type>0) THEN 531 IF(rayleigh_friction_type<99) THEN 517 532 phi=f_geopot(ind) 518 533 ue=f_ue(ind) … … 524 539 ENDDO 525 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 526 559 END IF 527 560 END DO … … 529 562 CALL trace_end("dissip") 530 563 564 CALL write_dissip_tendencies 531 565 !$OMP BARRIER 532 566 533 567 CONTAINS 568 534 569 SUBROUTINE relax(shift_t, shift_u) 535 570 USE dcmip_initial_conditions_test_1_2_3 … … 559 594 END SUBROUTINE relax 560 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 561 612 END SUBROUTINE dissip 562 613 -
codes/icosagcm/devel/src/icosa_init.f90
r686 r714 21 21 USE etat0_mod 22 22 USE diagflux_mod 23 USE profiling_mod 23 24 IMPLICIT NONE 24 25 26 CALL init_profiling 25 27 CALL init_mpipara 26 28 CALL trace_off … … 56 58 57 59 CALL init_diagflux 60 CALL zero_du_phys 58 61 CALL timeloop 59 62 CALL switch_omp_no_distrib_level -
codes/icosagcm/devel/src/output/output_field.f90
r533 r714 1 1 MODULE output_field_mod 2 USE genmod 2 USE genmod 3 USE xios_mod 4 USE profiling_mod 5 IMPLICIT NONE 6 SAVE 3 7 PRIVATE 8 4 9 LOGICAL,SAVE :: xios_output 5 10 !$OMP THREADPRIVATE(xios_output) 6 11 LOGICAL,SAVE :: enable_io 7 12 !$OMP THREADPRIVATE(enable_io) 13 14 INTEGER :: id_output 8 15 9 16 PUBLIC enable_io,xios_output,output_field_init,output_field,output_field_finalize … … 13 20 SUBROUTINE output_field_init 14 21 USE getin_mod 15 USE xios_mod16 USE write_field_mod17 22 IMPLICIT NONE 23 24 CALL register_id('output',id_output) 18 25 19 26 enable_io=.TRUE. … … 30 37 CALL xios_init_write_field 31 38 ENDIF 32 33 34 39 END SUBROUTINE output_field_init 35 40 36 41 SUBROUTINE output_field(name_in,field) 37 42 USE field_mod 38 USE xios_mod39 43 USE write_field_mod 40 44 IMPLICIT NONE 41 CHARACTER(LEN=*),INTENT(IN) :: name_in 42 TYPE(t_field),POINTER :: field(:) 43 44 IF (xios_output) THEN 45 CALL xios_write_field(name_in,field) 46 ELSE 45 CHARACTER(LEN=*),INTENT(IN) :: name_in 46 TYPE(t_field),POINTER :: field(:) 47 48 CALL enter_profile(id_output) 49 IF (xios_output) THEN 50 CALL xios_write_field(name_in,field) 51 ELSE 47 52 CALL writeField(name_in,field) 48 ENDIF 53 ENDIF 54 CALL exit_profile(id_output) 49 55 50 56 END SUBROUTINE output_field 51 57 52 58 SUBROUTINE output_field_finalize 53 USE ioipsl 54 USE xios_mod 55 IMPLICIT NONE 59 USE ioipsl 60 IMPLICIT NONE 56 61 57 62 IF (xios_output) THEN 58 CALL xios_write_field_finalize63 CALL xios_write_field_finalize 59 64 ENDIF 60 65 -
codes/icosagcm/devel/src/parallel/mpipara.F90
r533 r714 12 12 INTEGER,SAVE :: mpi_master 13 13 14 14 INTEGER,SAVE :: id_mpi ! id for profiling 15 15 16 INTERFACE allocate_mpi_buffer 16 17 MODULE PROCEDURE allocate_mpi_buffer_r2, allocate_mpi_buffer_r3,allocate_mpi_buffer_r4 -
codes/icosagcm/devel/src/parallel/transfert_mpi.f90
r603 r714 2 2 USE genmod 3 3 USE field_mod 4 IMPLICIT NONE 4 5 5 6 TYPE array … … 85 86 END INTERFACE 86 87 87 88 89 88 CONTAINS 90 89 91 90 92 91 SUBROUTINE init_transfert 92 USE profiling_mod 93 93 USE domain_mod 94 94 USE dimensions … … 100 100 INTEGER :: ind,i,j 101 101 LOGICAL ::ok 102 102 103 CALL register_id('MPI', id_mpi) 104 103 105 CALL create_request(field_t,req_i1) 104 106 … … 1094 1096 1095 1097 SUBROUTINE send_message_mpi(field,message) 1098 USE profiling_mod 1096 1099 USE field_mod 1097 1100 USE domain_mod … … 1124 1127 1125 1128 ! CALL trace_start("send_message_mpi") 1129 1130 CALL enter_profile(id_mpi) 1126 1131 1127 1132 !$OMP BARRIER … … 1482 1487 !$OMP BARRIER 1483 1488 ! CALL trace_end("send_message_mpi") 1489 1490 CALL exit_profile(id_mpi) 1484 1491 1485 1492 END SUBROUTINE send_message_mpi … … 1499 1506 1500 1507 SUBROUTINE wait_message_mpi(message) 1508 USE profiling_mod 1501 1509 USE field_mod 1502 1510 USE domain_mod … … 1527 1535 message%open=.FALSE. 1528 1536 IF (.NOT. message%pending) RETURN 1537 1538 CALL enter_profile(id_mpi) 1529 1539 1530 1540 ! CALL trace_start("wait_message_mpi") … … 1667 1677 !$OMP BARRIER 1668 1678 1679 CALL exit_profile(id_mpi) 1680 1669 1681 END SUBROUTINE wait_message_mpi 1670 1682 -
codes/icosagcm/devel/src/physics/physics.f90
r584 r714 1 1 MODULE physics_mod 2 2 USE icosa 3 3 USE field_mod 4 4 USE physics_interface_mod 5 USE omp_para 6 IMPLICIT NONE 5 7 PRIVATE 6 8 … … 14 16 TYPE(t_field),POINTER,SAVE :: f_p(:), f_pk(:) 15 17 TYPE(t_field),POINTER,SAVE :: f_temp(:) 18 TYPE(t_field),POINTER,SAVE :: f_du_phys(:) 16 19 17 20 CHARACTER(LEN=255),SAVE :: physics_type 18 21 !$OMP THREADPRIVATE(physics_type) 19 22 20 PUBLIC :: physics, init_physics 23 PUBLIC :: physics, init_physics, zero_du_phys 21 24 22 25 CONTAINS … … 25 28 USE mpipara 26 29 USE etat0_mod 27 USE icosa28 USE physics_interface_mod29 30 USE physics_dcmip_mod, ONLY : init_physics_dcmip=>init_physics 30 31 USE physics_dcmip2016_mod, ONLY : init_physics_dcmip2016=>init_physics … … 32 33 USE physics_lmdz_generic_mod, ONLY : init_physics_lmdz_generic=>init_physics 33 34 USE physics_external_mod, ONLY : init_physics_external=>init_physics 34 IMPLICIT NONE35 35 36 36 physics_inout%dt_phys = dt*itau_physics … … 84 84 END SELECT 85 85 86 CALL allocate_field(f_du_phys,field_u,type_real,llm, name='du_phys') 87 86 88 IF(is_mpi_root) PRINT *, 'phys_type = ',phys_type 87 89 END SUBROUTINE init_physics 88 90 91 SUBROUTINE zero_du_phys() 92 REAL(rstd), DIMENSION(:,:), POINTER :: du 93 INTEGER :: ind 94 DO ind=1,ndomain 95 IF (.NOT. assigned_domain(ind)) CYCLE 96 CALL swap_dimensions(ind) 97 CALL swap_geometry(ind) 98 du=f_du_phys(ind) 99 du(:,ll_begin:ll_end) = 0. 100 END DO 101 END SUBROUTINE zero_du_phys 102 103 SUBROUTINE add_du_phys(coef, f_u) 104 REAL(rstd), INTENT(IN) :: coef ! -1 before physics, +1 after physics 105 TYPE(t_field),POINTER :: f_u(:) ! velocity field before/after call to physics 106 REAL(rstd), DIMENSION(:,:), POINTER :: u, du 107 INTEGER :: ind 108 DO ind=1,ndomain 109 IF (.NOT. assigned_domain(ind)) CYCLE 110 CALL swap_dimensions(ind) 111 CALL swap_geometry(ind) 112 du=f_du_phys(ind) 113 u=f_u(ind) 114 du(:,ll_begin:ll_end) = du(:,ll_begin:ll_end) + coef*u(:,ll_begin:ll_end) 115 END DO 116 END SUBROUTINE add_du_phys 117 89 118 SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q) 90 USE icosa91 USE physics_interface_mod92 119 USE physics_lmdz_generic_mod, ONLY : physics_lmdz_generic => physics 93 120 USE physics_external_mod, ONLY : physics_external => physics … … 96 123 USE etat0_heldsz_mod 97 124 USE etat0_venus_mod, ONLY : phys_venus => physics 98 IMPLICIT NONE99 125 INTEGER, INTENT(IN) :: it 100 126 TYPE(t_field),POINTER :: f_phis(:) … … 109 135 TYPE(t_physics_inout) :: args 110 136 111 IF(MOD(it,itau_physics)==0) THEN 112 137 IF(MOD(it,itau_physics)==0 .AND. phys_type/=phys_none) THEN 138 139 ! as a result of the the two calls to add_du_phys, 140 ! du_phys increases by u(after physics) - u (before physics) 141 CALL add_du_phys(-1., f_ue) 142 113 143 SELECT CASE(phys_type) 114 CASE (phys_none)115 ! No physics, do nothing116 144 CASE(phys_HS94) 117 145 CALL held_suarez(f_ps,f_theta_rhodz,f_ue) … … 129 157 CALL transfert_request(f_ue,req_e0_vect) 130 158 CALL transfert_request(f_q,req_i0) 159 160 CALL add_du_phys(1., f_ue) 131 161 END IF 132 162 133 163 IF (mod(it,itau_out)==0 ) THEN 164 CALL write_physics_tendencies 165 CALL zero_du_phys 134 166 SELECT CASE(phys_type) 135 167 CASE (phys_DCMIP) … … 142 174 END SUBROUTINE physics 143 175 176 SUBROUTINE write_physics_tendencies 177 USE observable_mod, ONLY : f_buf_ulon, f_buf_ulat 178 USE wind_mod 179 USE output_field_mod 180 CALL transfert_request(f_du_phys,req_e1_vect) 181 CALL un2ulonlat(f_du_phys, f_buf_ulon, f_buf_ulat, (1./(dt*itau_out))) 182 CALL output_field("dulon_phys",f_buf_ulon) 183 CALL output_field("dulat_phys",f_buf_ulat) 184 END SUBROUTINE write_physics_tendencies 185 144 186 SUBROUTINE physics_column(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q) 145 USE icosa146 USE physics_interface_mod147 187 USE physics_dcmip_mod, ONLY : full_physics_dcmip => full_physics 148 188 USE physics_dcmip2016_mod, ONLY : full_physics_dcmip2016 => full_physics 149 189 USE theta2theta_rhodz_mod 150 190 USE mpipara 151 USE omp_para152 191 USE checksum_mod 153 IMPLICIT NONE154 192 TYPE(t_field),POINTER :: f_phis(:) 155 193 TYPE(t_field),POINTER :: f_ps(:) … … 229 267 230 268 SUBROUTINE pack_physics(info, phis, ps, temp, ue, q, p, pk, ulon, ulat ) 231 USE icosa232 269 USE wind_mod 233 270 USE pression_mod 234 271 USE theta2theta_rhodz_mod 235 USE physics_interface_mod236 272 USE exner_mod 237 USE omp_para238 IMPLICIT NONE239 273 TYPE(t_pack_info) :: info 240 274 REAL(rstd) :: phis(iim*jjm) … … 272 306 273 307 SUBROUTINE unpack_physics(info, ps,temp, q, dulon, dulat) 274 USE icosa275 USE physics_interface_mod276 308 USE theta2theta_rhodz_mod 277 USE omp_para278 IMPLICIT NONE279 309 TYPE(t_pack_info) :: info 280 310 REAL(rstd) :: ps(iim*jjm) … … 303 333 304 334 SUBROUTINE compute_update_velocity(dulon, dulat, ue) 305 USE icosa306 USE physics_interface_mod307 335 USE wind_mod 308 USE omp_para309 IMPLICIT NONE310 336 REAL(rstd) :: dulon(iim*jjm,llm) 311 337 REAL(rstd) :: dulat(iim*jjm,llm) -
codes/icosagcm/devel/src/time/timeloop_gcm.f90
r609 r714 1 1 MODULE timeloop_gcm_mod 2 USE profiling_mod 2 3 USE icosa 3 4 USE disvert_mod … … 13 14 TYPE(t_message),SAVE :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q0, req_W0, req_geopot0 14 15 LOGICAL, SAVE :: positive_theta 15 16 INTEGER :: itau_prof, id_timeloop, id_dyn, id_phys, id_dissip, id_adv, id_diags 16 17 PUBLIC :: init_timeloop, timeloop 17 18 … … 32 33 CHARACTER(len=255) :: def 33 34 35 CALL register_id('timeloop',id_timeloop) 36 CALL register_id('dyn',id_dyn) 37 CALL register_id('dissip',id_dissip) 38 CALL register_id('phys',id_phys) 39 CALL register_id('adv',id_adv) 40 CALL register_id('diags',id_diags) 41 34 42 CALL init_caldyn 35 43 36 44 ! IF (xios_output) itau_out=1 37 45 IF (.NOT. enable_io) itau_out=HUGE(itau_out) 46 47 itau_prof=100 48 CALL getin('itau_profiling',itau_prof) 38 49 39 50 positive_theta=.FALSE. … … 231 242 232 243 !$OMP MASTER 233 CALL SYSTEM_CLOCK(start_clock) 234 CALL SYSTEM_CLOCK(count_rate=rate_clock) 244 CALL SYSTEM_CLOCK(start_clock, rate_clock) 235 245 !$OMP END MASTER 236 246 237 247 DO it=itau0+1,itau0+itaumax 238 239 248 IF (is_master) CALL print_iteration(it, itau0, itaumax, start_clock, rate_clock) 249 250 CALL enter_profile(id_timeloop) 240 251 241 252 IF (xios_output) THEN … … 267 278 CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q) 268 279 280 CALL enter_profile(id_dyn) 269 281 SELECT CASE(scheme_family) 270 282 CASE(explicit) … … 273 285 CALL HEVI_scheme(it, fluxt_zero) 274 286 END SELECT 275 287 CALL exit_profile(id_dyn) 288 289 CALL enter_profile(id_dissip) 276 290 IF (MOD(it,itau_dissip)==0) THEN 277 291 … … 288 302 ENDIF 289 303 304 CALL enter_profile(id_diags) 290 305 CALL check_conserve_detailed(it, AAM_dyn, & 291 306 f_ps,f_dps,f_u,f_theta_rhodz,f_phis) 307 CALL exit_profile(id_diags) 292 308 293 309 CALL dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_u, f_dtheta_rhodz,f_du) … … 299 315 END IF 300 316 317 CALL enter_profile(id_diags) 301 318 CALL check_conserve_detailed(it, AAM_dissip, & 302 319 f_ps,f_dps,f_u,f_theta_rhodz,f_phis) 320 CALL exit_profile(id_diags) 303 321 END IF 322 CALL exit_profile(id_dissip) 304 323 324 CALL enter_profile(id_adv) 305 325 IF(MOD(it,itau_adv)==0) THEN 306 326 CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz, & ! update q and rhodz after RK step … … 322 342 IF(positive_theta) CALL copy_q_to_theta(f_theta_rhodz,f_rhodz,f_q) 323 343 END IF 344 CALL exit_profile(id_adv) 324 345 346 CALL enter_profile(id_diags) 325 347 IF (MOD(it,itau_physics)==0) THEN 326 348 CALL check_conserve_detailed(it, AAM_dyn, & 327 349 f_ps,f_dps,f_u,f_theta_rhodz,f_phis) 350 CALL enter_profile(id_phys) 328 351 CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q) 352 CALL exit_profile(id_phys) 329 353 CALL check_conserve_detailed(it, AAM_phys, & 330 354 f_ps,f_dps,f_u,f_theta_rhodz,f_phis) … … 345 369 CALL write_output_fields_basic(.FALSE.,f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q) 346 370 ENDIF 347 371 CALL exit_profile(id_diags) 372 373 CALL exit_profile(id_timeloop) 348 374 END DO 349 375 … … 382 408 ' -- Completion in (min) : ', INT((total-elapsed)/60.) 383 409 END IF 410 IF(MOD(it,itau_prof)==0) CALL print_profile 411 384 412 END SUBROUTINE print_iteration 385 413
Note: See TracChangeset
for help on using the changeset viewer.