!############################# Change Log ################################## ! 5.2.0 ! !MB: new file for the Runge-Kutta dynamical core ! (after Wicker, Skamarock, 2002, MWR) ! Saulo Freitas (INPE), Michael Baldauf (DWD) ! !############################# Change Log ################################## module ModTimestep_RK !MB: for testing only logical :: flag_mb_adv_test=.false. character(len=*), parameter :: h="**(rtimh_rk)**" contains subroutine timestep_rk(OneGrid) use ModMessageSet, only: & PostRecvSendMsgs, & WaitRecvMsgs use ModAcoust, only: & acoustic_new, & init_div_damping_coeff, & deallocate_alpha_div, & apply_div_damping ! use ModGrid, only: & Grid use mem_basic, only: & basic_g ! INTENT(INOUT) use node_mod, only: & mzp, mxp, myp, & ! INTENT(IN) ia, iz, ja, jz, & ! INTENT(IN) i0, j0, & ! INTENT(IN) izu, jzv, & ! INTENT(IN) mynum, & ! INTENT(IN) ibcon, & ! INTENT(IN) nmachs ! INTENT(IN) use mem_cuparm, only: & NNQPARM, & ! INTENT(IN) IF_CUINV ! INTENT(IN) use mem_varinit, only: & NUD_TYPE ! INTENT(IN) use mem_turb, only: & IF_URBAN_CANOPY, & ! INTENT(IN) ihorgrad ! INTENT(IN) use mem_oda, only: & if_oda ! INTENT(IN) use micphys, only: & mcphys_type, &! INTENT(IN) level ! INTENT(IN) use mem_grid, only: & ngrids, & ! INTENT(IN) ngrid, & ! INTENT(IN) npatch, & ! INTENT(IN) time, & ! INTENT(IN) dts, & dtlong, & ! INTENT(IN) dtlongn, & ! INTENT(IN) ideltat, & nnacoust, & iyear1, & ! INTENT(IN) imonth1, & ! INTENT(IN) idate1, & ! INTENT(IN) grid_g, & ! INTENT(INOUT) nxtnest, & ! INTENT(IN) if_adap, & ! INTENT(IN) dtlt, & ! INTENT(IN) istp, & ! INTENT(IN) jdim, & ! INTENT(IN) nzp, & ! INTENT(IN) f_thermo_e, & ! INTENT(IN) f_thermo_w, & ! INTENT(IN) f_thermo_s, & ! INTENT(IN) f_thermo_n, & ! INTENT(IN) zt, & zm, & dzt, & nzpmax, & itime1, & vveldamp ! INTENT(IN) use shcu_vars_const, only: & ! For Shallow Cumulus Paramet. NNSHCU ! INTENT(IN) use mem_scalar, only: & ! For SiB scalar_g ! INTENT(IN) use mem_leaf, only: & ! For SiB ISFCL ! INTENT(IN) ! TEB_SPM use teb_spm_start, only: & TEB_SPM ! INTENT(IN) use mem_emiss, only: & ichemi, & ! INTENT(IN) isource ! INTENT(IN) !ALF- Necessary in new advection scheme use advect_kit, only : & calc_advec ! Subroutine ! For specific optimization depending the type of machine use machine_arq, only: & machine ! INTENT(IN) use rconstants, only: & g, & ! (IN) cp, & ! (IN) cpor, & ! (IN) p00, & ! (IN) rgas, & ! (IN) pi180 ! (IN) use ccatt_start, only: & ccatt ! (IN) use mem_tend, only: & tend use mem_chem1, only: & nvert_src=>chem1_src_z_dim_g, & ! (IN) chem1_g, & ! (INOUT) nsrc, & ! (IN) chem1_src_g, & ! %sc_src(INOUT) chemistry, & ! (IN) split_method, & ! (IN) n_dyn_chem, & ntimes_src use mem_aer1, only: & aerosol, &! (IN) aer1_g, &! %sc_src(INOUT) aer2_g, &! %sc_src(INOUT) aer_nvert_src=>aer1_src_z_dim_g ! (IN) use mem_plume_chem1, only: plume_mean_g & ! %flam_frac(IN), %fire_size(IN) ,plume_fre_g use mem_stilt, only: & iexev, & ! (IN) imassflx, & ! (IN) stilt_g ! %dnp (IN) use mem_radiate, only: radiate_g use chem_sources, only : & alloc_emiss_cycle, & ! Subroutine init_actual_time_index, & ! Subroutine emiss_cycle, & ! (INOUT) emiss_cycle_alloc, & srcmapfn ! (IN) use ChemSourcesDriver, only: sources_driver ! Subroutine use ChemDryDepDriver , only: drydep_driver ! Subroutine use module_chemistry_driver, only: chemistry_driver ! Subroutine use radiation, only: radiate ! Subroutine use ModTimeStamp, only: SynchronizedTimeStamp, TimeStamp use cuparm_grell3, only: cuparm_grell3_catt & ! subroutine ,g3d_g use digitalFilter, only: & applyDigitalFilter, & ! subroutine fileNameDF, & ! intent(inout) - file control dfVars, & applyDF USE monotonic_adv, only: & advmnt_driver, & ! subroutine advmnt USE DriverMatrix, ONLY: MatrixDriver !Matrix Aerosol Model USE rams_microphysics_2M, only: micro_2M_rams60,negadj1_2M_rams60 use rrtm_driv, only: rrtm_driver USE mem_radiate, ONLY: & ilwrtyp, iswrtyp use CUPARM_GRELL3, only: g3d_g use mem_turb, only: & turb_g use ModComm, only: commHaloAcou use wind_farm, only: wind_farm_driver,windfarm implicit none type(Grid), pointer :: OneGrid ! execution time instrumentation include "tsNames.h" include "i8.h" LOGICAL, parameter :: flag_Coriolis_in_every_RK_step = .FALSE. integer :: l_rk real :: rk_beta(3) integer :: rk_nmbr_small_timesteps(3) integer, parameter :: &!rk_order = 2 ! for testing purposes (but works only with upwind1 and upwind3 advection) rk_order = 3 ! for Wicker, Skamarock (2002) MWR-scheme integer :: i,j,k integer :: n character(len=2) :: crk logical :: singleProcRun singleProcRun = nmachs == 1 !MB: only for testing !integer :: nmbr_gpts !real :: pm,tm ! +------------------------------------------------------------------+ ! | Timestep driver for the Runge-Kutta non-hydrostatic time-split | ! | model. | ! +------------------------------------------------------------------+ ! Zero out all tendency arrays. !-------------------------------- call TEND0() if (iexev == 2) & call exevolve(mzp,mxp,myp,ngrid,ia,iz,ja,jz,izu,jzv,jdim,mynum,dtlt,'ADV') ! Send boundaries to adjoining nodes !------------------------------------------- if (nmachs > 1) then call PostRecvSendMsgs(OneGrid%SelectedGhostZoneSend, OneGrid%SelectedGhostZoneRecv) endif ! Analysis nudging and boundary condition !------------------------------------------ if (NUD_TYPE>0) call DATASSIM() ! Get the overlap region between parallel nodes !--------------------------------------------------- if (nmachs > 1) then call WaitRecvMsgs(OneGrid%SelectedGhostZoneSend, OneGrid%SelectedGhostZoneRecv) endif if (iexev == 2) & call exevolve(mzp,mxp,myp,ngrid,ia,iz,ja,jz,izu,jzv,jdim,mynum,dtlt,'THA') if (iexev == 2) & call exevolve(mzp,mxp,myp,ngrid,ia,iz,ja,jz,izu,jzv,jdim,mynum,dtlt,'THS') !------------------------------------------------------------------------------ ! init preparations for Runge-Kutta -loop !------------------------------------------------------------------------------ if ( rk_order == 2 ) then ! Wicker, Skamarock (1998)-RK-scheme rk_beta(1) = 1.0 / 2.0 ! = beta(2,1) of Butcher tableau rk_beta(2) = 1.0 ! = beta(3,2) of Butcher tableau if ( MOD( nnacoust(ngrid), 2) /= 0 ) then call fatal_error("ERROR in timestep_rk: nnacoust(ngrid) must be an integer multiple of 2") end if rk_nmbr_small_timesteps(1) = nnacoust(ngrid) / 2 rk_nmbr_small_timesteps(2) = nnacoust(ngrid) else if ( rk_order == 3 ) then ! Wicker, Skamarock (2002)-RK-scheme rk_beta(1) = 1.0 / 3.0 ! = beta(2,1) of Butcher tableau rk_beta(2) = 1.0 / 2.0 ! = beta(3,2) of Butcher tableau rk_beta(3) = 1.0 ! = beta(4,3) of Butcher tableau if ( MOD( nnacoust(ngrid), 6) /= 0 ) then call fatal_error("ERROR in timestep_rk: nnacoust(ngrid) must be an integer multiple of 6") end if !- no futuro, tente rk_nmbr_small_timesteps(1)=1 e dts=2* dtlt / nnacoust(ngrid) rk_nmbr_small_timesteps(1) = nnacoust(ngrid) / 3 rk_nmbr_small_timesteps(2) = nnacoust(ngrid) / 2 rk_nmbr_small_timesteps(3) = nnacoust(ngrid) else call fatal_error("ERROR in timestep_rk: false value for rk_order") end if dts = dtlt / nnacoust(ngrid) if ( apply_div_damping ) then ! calculation of divergence damping coeff: !MB: ATTENTION: dts different for different nests?!! if ( ideltat == 0 ) then ! it is sufficient to calculate alpha_div only once: if ( istp == 1 ) then call init_div_damping_coeff( dts ) end if else call init_div_damping_coeff( dts ) end if end if ! begin of Runge-Kutta loop !--------------------------------------------------- !MB>> !nmbr_gpts = mxp * myp * mzp !MB: only for testing!!! !write(*,"(A,I2,4F15.8)") "uc ", 0, minval( basic_g(ngrid)%uc ), maxval( basic_g(ngrid)%uc ), & ! sum( basic_g(ngrid)%uc )/nmbr_gpts, basic_g(ngrid)%uc(7,8,9) !write(*,"(A,I2,4F15.8)") "vc ", 0, minval( basic_g(ngrid)%vc ), maxval( basic_g(ngrid)%vc ), & ! sum( basic_g(ngrid)%vc )/nmbr_gpts, basic_g(ngrid)%vc(7,8,9) !write(*,"(A,I2,4F15.8)") "wc ", 0, minval( basic_g(ngrid)%wc ), maxval( basic_g(ngrid)%wc ), & ! sum( basic_g(ngrid)%wc )/nmbr_gpts, basic_g(ngrid)%wc(7,8,9) !write(*,"(A,I2,4F15.8)") "pc ", 0, minval( basic_g(ngrid)%pc ), maxval( basic_g(ngrid)%pc ), & ! sum( basic_g(ngrid)%pc )/nmbr_gpts, basic_g(ngrid)%pc(7,8,9) !write(*,"(A,I2,4F15.8)") "thc ", 0, minval( basic_g(ngrid)%thc ), maxval( basic_g(ngrid)%thc ), & ! sum( basic_g(ngrid)%thc )/nmbr_gpts, basic_g(ngrid)%thc(7,8,9) ! !write(*,"(A,I2,3F15.8)") "ut ", 0, minval( tend%ut ), maxval( tend%ut ), sum( tend%ut )/nmbr_gpts !write(*,"(A,I2,3F15.8)") "vt ", 0, minval( tend%vt ), maxval( tend%vt ), sum( tend%vt )/nmbr_gpts !write(*,"(A,I2,3F15.8)") "wt ", 0, minval( tend%wt ), maxval( tend%wt ), sum( tend%wt )/nmbr_gpts !write(*,"(A,I2,3F15.8)") "pt ", 0, minval( tend%pt ), maxval( tend%pt ), sum( tend%pt )/nmbr_gpts !write(*,"(A,I2,3F15.8)") "tht ", 0, minval( tend%tht ), maxval( tend%tht ), sum( tend%tht )/nmbr_gpts !!MB<< !MB: !if ( flag_mb_adv_test) then ! !! ascii output for horizontal scalar advection test: ! !do i=1, mxp ! ! write(*,*) "thc ", i, basic_g(ngrid)%thc(6,i,5) ! !end do ! !! ascii output for vertical scalar advection test: ! do k=1, mzp ! write(*,*) "thc ", k, zt(k), basic_g(ngrid)%thc(k,10,5) ! end do !end if ! Lateral velocity boundaries - radiative !------------------------------------------- call LATBND() do l_rk = 1, rk_order !write(crk,fmt='(I2.2)') l_rk !print*,"rk:",rk_beta(l_rk),rk_nmbr_small_timesteps(l_rk),& ! rk_nmbr_small_timesteps(l_rk)*dts,dtlt !initialize the tendencies with the physics tendencies tend%ut_rk (:) =tend%ut (:) tend%vt_rk (:) =tend%vt (:) tend%wt_rk (:) =tend%wt (:) tend%pt_rk (:) =tend%pt (:) !; call dumpVarAllLatLonk(tend%tht,'THT',515,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0) tend%tht_rk(:) =tend%tht(:) !; call dumpVarAllLatLonk(tend%tht_rk,'THT_RK',516,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0) !.. call dumpVarAllLatLonk3P(tend%ut_rk , 'UT' ,528,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0,h) !.. call dumpVarAllLatLonk3P(tend%vt_rk , 'VT' ,528,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0,h) !.. call dumpVarAllLatLonk3P(tend%wt_rk , 'WT' ,528,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0,h) !.. call dumpVarAllLatLonk3P(tend%pt_rk , 'PT' ,528,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0,h) !.. call dumpVarAllLatLonk3P(tend%tht_rk,'THT' ,528,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0,h) ! advection should give back tendencies ! ut_rk, vt_rk, wt_rk, pt_rk, tht_rk = physics tend + advection tendency ! Velocity advection !---------------------------------------- call advectc_rk('V',mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum,l_rk)!; call dumpVarAllLatLonk(tend%wt_rk , 'WT' ,530,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0) ! advection of pi and theta_il !---------------------------------------- !call dumpVarAllLatLonk(basic_g(ngrid)%thp,'aTHPe'//crk,1,mxp,1,myp,1,mzp,0.0,0.0) !ok na 1 !call dumpVarAllLatLonk(basic_g(ngrid)%thc,'aTHCe'//crk,1,mxp,1,myp,1,mzp,0.0,0.0) !ok na 1 call advectc_rk('THETAIL',mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum,l_rk)!; call dumpVarAllLatLonk(tend%wt_rk , 'WT' ,536,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0) call advectc_rk('PI' ,mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum,l_rk)!; call dumpVarAllLatLonk(tend%wt_rk , 'WT' ,537,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0) ! call dumpVarAllLatLonk(tend%tht_rk,'a1THT'//crk,1,mxp,1,myp,1,mzp,0.0,0.0) !furado na 1 ! if ( flag_Coriolis_in_every_RK_step ) then !if ( .not. flag_mb_adv_test) call CORLOS(mzp,mxp,myp,i0,j0,ia,iz,ja,jz,izu,jzv, tend%ut_rk, tend%vt_rk) !end if end if !call dumpVarAllLatLonk(tend%wt_rk , 'WT' ,553,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0) ! Buoyancy term for w equation !---------------------------------------- !if ( .not. flag_mb_adv_test) call BUOYANCY( tend%wt_rk ) !end if !call dumpVarAllLatLonk(tend%wt_rk , 'WT' ,559,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0) if ( l_rk > 1 ) then ! (not necessary in the first RK substep) basic_g(ngrid)%uc (:,:,:) = basic_g(ngrid)%up (:,:,:) basic_g(ngrid)%vc (:,:,:) = basic_g(ngrid)%vp (:,:,:) basic_g(ngrid)%wc (:,:,:) = basic_g(ngrid)%wp (:,:,:) basic_g(ngrid)%pc (:,:,:) = basic_g(ngrid)%pp (:,:,:) basic_g(ngrid)%thc(:,:,:) = basic_g(ngrid)%thp(:,:,:) end if !.. call dumpVarAllLatLonk3P(tend%ut_rk , 'UT' ,561,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0,h) !.. call dumpVarAllLatLonk3P(tend%vt_rk , 'VT' ,561,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0,h) !.. call dumpVarAllLatLonk3P(tend%wt_rk , 'WT' ,561,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0,h) !.. call dumpVarAllLatLonk3P(tend%pt_rk , 'PT' ,561,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0,h) !.. call dumpVarAllLatLonk3P(tend%tht_rk,'THT' ,561,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0,h) !- Acoustic small timesteps if ( .not. flag_mb_adv_test) then call acoustic_new(OneGrid, rk_nmbr_small_timesteps(l_rk),l_rk ) !call acoustic_new(OneGrid, rk_nmbr_small_timesteps(l_rk),l_rk ) else ! just for testing: ! without the acoustic subroutine, tendencies must be added "manually" ! update pp -> pc (similar for u, v, w, if needed.) call update_long_rk(int(mxp*myp*mzp,i8),dtlt,rk_beta(l_rk) & ,basic_g(ngrid)%pc,basic_g(ngrid)%pp & ,tend%pt_rk) endif !.. call dumpVarAllLatLonk3P(basic_g(ngrid)%uc, 'UC' ,586,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0,h) !.. call dumpVarAllLatLonk3P(basic_g(ngrid)%vc, 'VC' ,586,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0,h) !.. call dumpVarAllLatLonk3P(basic_g(ngrid)%wc, 'WC' ,586,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0,h) !.. call dumpVarAllLatLonk3P(basic_g(ngrid)%pc, 'PC' ,586,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0,h) !- update thp -> thc (theta_il is not contained in acoustic_new) if (.not. singleProcRun) then call commHaloAcou(tend%tht_rk,mxp,myp,mzp,myNum,'tht_rk') endif call update_long_rk(int(mxp*myp*mzp,i8),dtlt,rk_beta(l_rk) & ,basic_g(ngrid)%thc,basic_g(ngrid)%thp & ,tend%tht_rk) !call dumpVarAllLatLonk(basic_g(ngrid)%thc,'THC' ,602,l_rk,0,1,mxp,1,myp,1,mzp,0.0,0.0,h) ! call dumpVarAllLatLonk(basic_g(ngrid)%thc,'aTHCo'//crk,1,mxp,1,myp,1,mzp,0.0,0.0) !- determine theta (dry potential temp.) for the buoyancy term: call theta_thp_rk(mzp,mxp,myp,ia,iz,ja,jz,"get_theta") !-damping on vertical velocity to keep stability !MB: does this act on wc??? if(vveldamp == 1) call w_damping(mzp,mxp,myp,ia,iz,ja,jz,mynum) end do ! end of Runge-Kutta loop ! i.e. the fields ! basic_g%uc, ..%vc, ..%wc, ..%pc, ..%thc ! contain the fields at the new time level n+1 !--------------------------------------------------- ! !.. call dumpVarAllLatLonk3P(basic_g(ngrid)%uc, 'UC' ,605,0,0,1,mxp,1,myp,1,mzp,0.0,120.0,h) !.. call dumpVarAllLatLonk3P(basic_g(ngrid)%vc, 'VC' ,605,0,0,1,mxp,1,myp,1,mzp,0.0,120.0,h) !.. call dumpVarAllLatLonk3P(basic_g(ngrid)%wc, 'WC' ,605,0,0,1,mxp,1,myp,1,mzp,0.0,120.0,h) !.. call dumpVarAllLatLonk3P(basic_g(ngrid)%pc, 'PC' ,605,0,0,1,mxp,1,myp,1,mzp,0.0,120.0,h) !.. call dumpVarAllLatLonk3P(basic_g(ngrid)%thc,'THC' ,605,0,0,1,mxp,1,myp,1,mzp,0.0,120.0,h) if ( apply_div_damping ) then if ( ideltat /= 0 ) then call deallocate_alpha_div end if end if ! ! ! water species, tke and tracers advection !---------------------------------------- IF(advmnt == 1) THEN !- monotonic advection scheme CALL advmnt_driver('SCALAR',mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum) ELSEIF(advmnt == 0) THEN !- using the 2nd order forward upstream CALL ADVECTc ('SCALAR',mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum) ELSEIF(advmnt == 3) THEN !- using the WS advection CALL advectc_rk('SCALAR',mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum,l_rk) ENDIF ! Update scalars (water species, tke and tracers) !---------------------------------------- ! call PREDTR() !- copy current time into past time (u,v,w,pi,thetail) !---> thp must be changed to THC for microphysics/bc/theta update !---> pp must be changed to PC for microphysics !---> wp must be changed to WC for microphysics !---> up,vp must be changed to UC,VC for output basic_g(ngrid)%up (:,:,:) = basic_g(ngrid)%uc (:,:,:) basic_g(ngrid)%vp (:,:,:) = basic_g(ngrid)%vc (:,:,:) basic_g(ngrid)%wp (:,:,:) = basic_g(ngrid)%wc (:,:,:) basic_g(ngrid)%pp (:,:,:) = basic_g(ngrid)%pc (:,:,:) basic_g(ngrid)%thp(:,:,:) = basic_g(ngrid)%thc(:,:,:) !----> !----> ! Apply scalar b.c.'s (THP is changed here) !---------------------------------------- call TRSETS() !---> THC must be changed to THP to include microphysics/trsets changes !---> for the next timestep basic_g(ngrid)%thc(:,:,:) = basic_g(ngrid)%thp(:,:,:) !---> ! Lateral velocity boundaries - radiative !------------------------------------------- !srf call LATBND() ! Velocity/pressure boundary conditions !---------------------------------------- call VPSETS() !- call THERMO on the boundaries call thermo_boundary_driver((time+dtlongn(ngrid)), dtlong, & f_thermo_e, f_thermo_w, f_thermo_s, f_thermo_n, & nzp, mxp, myp, jdim) end subroutine timestep_rk end module ModTimestep_RK !************************************************************************* !-srf- temp routine only for testing - not being used ! subroutine adv_p_driver(m1,m2,m3,ifm,ia,iz,ja,jz,izu,jzv,jdim,mynum,edt,key) use mem_basic, only: basic_g use mem_grid, only: grid_g, itopo,dyncore_flag use mem_stilt, only: stilt_g use mem_tend, only: tend use mem_scratch, only: scratch use micphys, only: level !if(vapour_on) use therm_lib, only: vapour_on implicit none !----- Arguments -----------------------------------------------------------------------! character(len=*) , intent(in) :: key integer , intent(in) :: m1,m2,m3,ifm,ia,iz,ja,jz,izu,jzv,jdim,mynum real , intent(in) :: edt !----- Local variables -----------------------------------------------------------------! integer :: i,j,k !---------------------------------------------------------------------------------------! call adv_p(m1,m2,m3,ia,iz,ja,jz,izu,jzv,jdim,itopo & ,grid_g(ifm)%rtgu ,grid_g(ifm)%fmapui & ,grid_g(ifm)%rtgv ,grid_g(ifm)%fmapvi & ,grid_g(ifm)%f13t ,grid_g(ifm)%f23t & ,grid_g(ifm)%rtgt ,grid_g(ifm)%fmapt & ,grid_g(ifm)%dxt ,grid_g(ifm)%dyt & ,basic_g(ifm)%uc ,basic_g(ifm)%dn0u & ,basic_g(ifm)%vc ,basic_g(ifm)%dn0v & ,basic_g(ifm)%dn0 ,basic_g(ifm)%wc & ,basic_g(ifm)%pc ,tend%pt ) end subroutine adv_p_driver subroutine adv_p(m1,m2,m3,ia,iz,ja,jz,izu,jzv,jdim,itopo,rtgu,fmapui,rtgv,fmapvi,f13t & ,f23t,rtgt,fmapt,dxt,dyt,uc,dn0u,vc,dn0v,dn0,wc,pc,pt) use mem_grid , only : hw4 & ! intent(in) , dzt ! ! intent(in) implicit none !----- Arguments -----------------------------------------------------------------------! integer , intent(in) :: m1,m2,m3,ia,iz,ja,jz,izu,jzv,itopo,jdim real , intent(in) , dimension(m2,m3) :: rtgu,fmapui,rtgv,fmapvi,f13t,f23t real , intent(in) , dimension(m2,m3) :: rtgt,fmapt,dxt,dyt real , intent(in) , dimension(m1,m2,m3) :: uc,dn0u,vc,dn0v,dn0,wc,pc real , intent(inout), dimension(m1,m2,m3) :: pt !----- Local variables -----------------------------------------------------------------! integer :: i,j,k,im,jm real :: c1z,c1x,c1y real , dimension(m1,m2,m3) :: flxu,flxv,flxw !---------------------------------------------------------------------------------------! !----- Compute momentum fluxes flxu, flxv, flxw ----------------------------------------! do j = 1,m3 do i = 1,m2 do k = 1,m1 flxu(k,i,j) = uc(k,i,j) * dn0u(k,i,j) * rtgu(i,j) * fmapui(i,j) flxv(k,i,j) = vc(k,i,j) * dn0v(k,i,j) * rtgv(i,j) * fmapvi(i,j) enddo enddo enddo if(itopo == 0) then do j = 1,m3 do i = 1,m2 do k = 1,m1-1 flxw(k,i,j) = wc(k,i,j) * .5 * (dn0(k,i,j) + dn0(k+1,i,j)) end do end do end do else do j = 1,m3 jm = max(j-1,1) do i = 1,m2 im = max(i-1,1) do k = 1,m1-1 flxw(k,i,j) = wc(k,i,j) * .5 * (dn0(k,i,j) + dn0(k+1,i,j)) & + hw4(k) * ( ( flxu(k,i,j) + flxu(k+1,i,j) & + flxu(k,im,j) + flxu(k+1,im,j) ) * f13t(i,j) & + ( flxv(k,i,j) + flxv(k+1,i,j) & + flxv(k,i,jm) + flxv(k+1,i,jm) ) * f23t(i,j) ) end do end do end do end if !---------------------------------------------------------------------------------------! ! Compute advection contribution of zonal gradient to Exner function tendency. ! !---------------------------------------------------------------------------------------! do j = ja,jz do i = ia,izu c1x = 0.5 / rtgt(i,j) * fmapt(i,j) * dxt(i,j) do k = 2,m1-1 pt(k,i,j) = pt(k,i,j) & - c1x / dn0(k,i,j) & * ( flxu(k,i,j) * (pc(k,i,j) + pc(k,i+1,j)) & - flxu(k,i-1,j) * (pc(k,i,j) + pc(k,i-1,j)) & - (flxu(k,i,j) - flxu(k,i-1,j)) * 2.* pc(k,i,j) ) end do end do end do !---------------------------------------------------------------------------------------! ! Compute advection contribution of meridional gradient to Exner function tendency. ! !---------------------------------------------------------------------------------------! do j=ja,jzv do i=ia,iz c1y = 0.5 / rtgt(i,j) * fmapt(i,j) * dyt(i,j) do k=2,m1-1 pt(k,i,j) = pt(k,i,j) & - c1y /dn0(k,i,j) & * ( flxv(k,i,j) * (pc(k,i,j)+pc(k,i,j+jdim)) & - flxv(k,i,j-jdim) * (pc(k,i,j)+pc(k,i,j-jdim)) & - (flxv(k,i,j)-flxv(k,i,j-jdim)) * 2.* pc(k,i,j) ) end do end do end do !---------------------------------------------------------------------------------------! ! Compute advection contribution of vertical gradient to Exner function tendency. ! !---------------------------------------------------------------------------------------! do j=ja,jz do i=ia,iz c1z = 0.5 / rtgt(i,j) do k=2,m1-1 pt(k,i,j) = pt(k,i,j) & - c1z * dzt(k) /dn0(k,i,j) & * ( flxw(k,i,j) * (pc(k,i,j)+pc(k+1,i,j)) & - flxw(k-1,i,j) * (pc(k,i,j)+pc(k-1,i,j)) & - (flxw(k,i,j)-flxw(k-1,i,j)) * 2. * pc(k,i,j) ) end do end do end do return end subroutine adv_p !==========================================================================================!