!### macro's ##################################################### ! #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if ! #include "tm5.inc" ! !################################################################# MODULE ADVECT use GO, only : gol, goPr, goErr implicit none private public :: Advect_Init, Advect_Done public :: dynam0 ! --- const -------------------------------------- character(len=*), parameter :: mname = 'advect' CONTAINS ! ================================================================ subroutine Advect_Init( status ) use dims , only : nregions use Meteo, only : Set use Meteodata, only : sp1_dat, sp2_dat, sp_dat #ifdef with_prism use MeteoData, only : tsp_dat #endif use MeteoData, only : mfu_dat, mfv_dat, mfw_dat use MeteoData, only : pu_dat, pv_dat, pw_dat ! --- in/out -------------------------------- integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/Advect_Init' ! --- local -------------------------------- integer :: n ! --- begin -------------------------------- ! loop over all (zoom) regions: do n = 1, nregions ! surface pressures and mass fluxes for advection call Set( sp_dat(n), status, used=.true. ) call Set( sp1_dat(n), status, used=.true. ) call Set( sp2_dat(n), status, used=.true. ) #ifdef with_prism call Set( tsp_dat(n), status, used=.true. ) #endif call Set( mfu_dat(n), status, used=.true. ) call Set( mfv_dat(n), status, used=.true. ) call Set( mfw_dat(n), status, used=.true. ) call Set( pu_dat(n), status, used=.true. ) call Set( pv_dat(n), status, used=.true. ) call Set( pw_dat(n), status, used=.true. ) end do ! regions ! ok status = 0 end subroutine Advect_Init ! *** subroutine Advect_Done( status ) ! --- in/out -------------------------------- integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/Advect_Done' ! --- begin -------------------------------- ! nothing to be done ! ok status = 0 end subroutine Advect_Done ! ================================================================ !----------------------------------------------------------------------- ! !**** dynam0 - calculates vertical massfluxes v 9.1 ! ! programmed by mh mpi HH 1-oct-1991 ! ! purpose ! ------- ! This sr is CALLed each time after new massfluxes are read in. ! Calculate the new air-masses in each grid box from the ! surface pressure. ! Calculate the amount of air moved in each advection substep, ! also in the vertical direction. ! ! interface ! --------- ! CALL dynam0 ! ! method ! ------ ! integrate air conservation equation in the vertical direction ! in order to obtain the vertical massfluxes. ! ! externals ! --------- ! none ! ! reference ! --------- ! see manual ! ! modified by mk (1999) for zoom version. ! 25 Jan 2012 - P. Le Sager - adapted for MPI lat-lon domain decomposition ! !----------------------------------------------------------------------- SUBROUTINE DYNAM0(region, status) use dims, only : im, jm, lm, ndyn, tref, zero, bt, xcyc use global_data, only : wind_dat use meteodata , only : pu_dat, pv_dat use tm5_distgrid, only : dgrid, Get_DistGrid, update_halo ! in/output integer,intent(in) :: region integer,intent(out) :: status ! local real, dimension(:,:,:), pointer :: pu,pv,am,bm,cm real, dimension(:,:), allocatable :: pit real, dimension(:,:,:), allocatable :: sd real, dimension(:,:,:), allocatable :: conv_adv real :: dtu, dtv, dtw integer :: i, j, l, lmr, i1, i2, j1, j2 character(len=*), parameter :: rname = mname//'/dynam0' ! start call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2 , J_STRT=j1, J_STOP=j2 ) allocate( pit(i1:i2, j1:j2 ) ) allocate( sd(i1:i2, j1:j2, lm(region)-1) ) allocate( conv_adv(i1:i2, j1:j2, lm(region) ) ) ! leave if fluxes are not in use: if ( .not. pu_dat(region)%used ) return if ( .not. pv_dat(region)%used ) return ! set pointers: pu => pu_dat(region)%data pv => pv_dat(region)%data am => wind_dat(region)%am bm => wind_dat(region)%bm cm => wind_dat(region)%cm !---- length of advection substeps (in seconds) in the three directions dtu=ndyn/(2.*tref(region)) ! again removed the revert! cmk dtv=ndyn/(2.*tref(region)) dtw=ndyn/(2.*tref(region)) ! number of levels in this region: lmr = lm(region) ! init to zero: am = zero bm = zero cm = zero ! update pu and pv before using them CALL UPDATE_HALO( dgrid(region), pu_dat(region)%data, pu_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) CALL UPDATE_HALO( dgrid(region), pv_dat(region)%data, pv_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) ! ! compute conv_adv, the horizontal mass convergence ! the arrays pu and pv contain the horizontal air mass fluxes crossing ! the grid box boundaries. pu(i,j,l) is the eastward mass flux in kg/sec ! at the eastern edge of box i,j,l; pv(i,j,l) is the northward ! mass flux at the southern edge of box i,j,l ! do l=1,lmr do j=j1, j2 do i=i1, i2 conv_adv(i,j,l)=pu(i-1,j,l)-pu(i,j,l)+pv(i,j,l)-pv(i,j+1,l) end do end do end do ! ! compute pit, the vertically integrated conv_adv ! do j=j1, j2 do i=i1, i2 pit(i,j)=conv_adv(i,j,1) do l=2,lmr pit(i,j)=pit(i,j)+conv_adv(i,j,l) end do ! ! compute the vertical massflux on the box boundaries (in kg/s) !mk in the hybrid system, the tendency in ps is given by bt. Bt is !mk already normalized. Now distribute pit such that the new mass !mk distribution after advection exactly matches the !mk exactly coincides with the distribution that is obtained !mk from the new surface pressure. ! sd(i,j,lmr-1)=conv_adv(i,j,lmr) - (bt(lmr)-bt(lmr+1))*pit(i,j) do l=lmr-2,1,-1 sd(i,j,l)=sd(i,j,l+1)+conv_adv(i,j,l+1)-(bt(l+1)-bt(l+2))*pit(i,j) end do end do end do ! ! compute amount of air moved each advection substep, am, bm or cm (kg) ! ! am(i,j,l) is the amount of air moved eastward at the eastern ! boundary of grid box i,j,l ! compute cm ! cm(i,j,l) is the amount of air moved in upward direction at the ! upper boundary of grid box i,j,l ! compute bm ! bm(i,j,l) is the amount of air moved northward at the southern ! boundary of grid box i,j,l ! do l= 1, lmr do j= j1, j2 do i= i1-1, i2 am(i,j,l)=dtu*pu(i,j,l) end do end do end do do l= 1, lmr do j= j1, j2+1 do i= i1, i2 bm(i,j,l)=dtv*pv(i,j,l) end do end do end do do l= 1, lmr-1 do j= j1, j2 do i= i1, i2 cm(i,j,l)=-dtw*sd(i,j,l) end do end do end do CALL UPDATE_HALO( dgrid(region), am, wind_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) nullify(pu) nullify(pv) nullify(am) nullify(bm) nullify(cm) deallocate( pit, sd, conv_adv) END SUBROUTINE DYNAM0 end module advect