!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| module surface_hgt !BOP ! !MODULE: surface_hgt ! ! !DESCRIPTION: ! Contains routine for computing change in surface height. ! ! !REVISION HISTORY: ! SVN:$Id: surface_hgt.F90 26603 2011-01-28 23:09:02Z njn01 $ ! !USES: use kinds_mod, only: int_kind, r8 use blocks, only: nx_block, ny_block, block, get_block ! use distribution, only: use domain_size use domain, only: nblocks_clinic, blocks_clinic use constants, only: grav, c0 use prognostic, only: max_blocks_clinic, PSURF, GRADPX, GRADPY, newtime, & curtime, oldtime use forcing_fields, only: FW, FW_OLD use grid, only: sfc_layer_type, sfc_layer_varthick, sfc_layer_rigid, & sfc_layer_oldfree, CALCU, tgrid_to_ugrid use time_management, only: mix_pass, dtp use tavg, only: define_tavg_field, accumulate_tavg_field, accumulate_tavg_now use movie, only: define_movie_field, movie_requested, update_movie_field implicit none private save ! !PUBLIC MEMBER FUNCTIONS: public :: init_surface_hgt, & dhdt !EOP !BOC !----------------------------------------------------------------------- ! ! tavg ids for surface height diagnostics ! !----------------------------------------------------------------------- integer (int_kind) :: & tavg_SSH, &! tavg id for sea surface height tavg_SSH2, &! tavg id for sea surface height squared (formerly H2) tavg_H3 ! tavg id for (Dx(SSH))**2 + (Dy(SSH))**2 !----------------------------------------------------------------------- ! ! movie ids for surface height diagnostics ! !----------------------------------------------------------------------- integer (int_kind) :: & movie_SSH ! movie id for sea surface height !EOC !*********************************************************************** contains !*********************************************************************** !BOP ! !IROUTINE: init_surface_hgt ! !INTERFACE: subroutine init_surface_hgt ! !DESCRIPTION: ! Initialized quantities related to SSH, namely some tavg diagnostic ! ids. ! ! !REVISION HISTORY: ! same as module !EOP !BOC !----------------------------------------------------------------------- ! ! define tavg diagnostic fields ! !----------------------------------------------------------------------- call define_tavg_field(tavg_SSH,'SSH',2, & long_name='Sea Surface Height', & units='centimeter', grid_loc='2110', & coordinates='TLONG TLAT time') ! formerly H2 call define_tavg_field(tavg_SSH2,'SSH2',2, & long_name='SSH**2', & units='cm^2', grid_loc='2110', & coordinates='TLONG TLAT time') call define_tavg_field(tavg_H3,'H3',2, & long_name='(Dx(SSH))**2 + (Dy(SSH))**2', & units='----', grid_loc='2110', & coordinates='TLONG TLAT time') !----------------------------------------------------------------------- ! ! define movie diagnostic fields ! !----------------------------------------------------------------------- call define_movie_field(movie_SSH,'SSH',0, & long_name='Sea Surface Height', & units='cm', grid_loc='2110') !----------------------------------------------------------------------- !EOC end subroutine init_surface_hgt !*********************************************************************** !BOP ! !IROUTINE: dhdt ! !INTERFACE: subroutine dhdt(DH,DHU) ! !DESCRIPTION: ! This routine calculates the change in surface height ! $d\eta/dt$ from surface pressure and including freshwater ! flux for a variable-thickness surface layer ! $(d\eta/dt - F_W)$. ! ! The surface pressure contribution at T points is ! \begin{equation} ! {{d\eta}\over{dt}} = {{(p^n - p^{n-1})}\over{g\Delta t}}. ! \end{equation} ! At U points, the change in surface height is the area-weighted ! average of its value at surrounding T points. ! ! In the advection routines, the vertical velocity $w = w_t,w_u$ in ! T,U columns is determined by integrating the continuity equation ! $L(1) = 0$ from the surface down to level k. Depending on surface ! layer type, the integration starts with: ! \begin{equation}\begin{array}{llll} ! w &=& 0 &{\rm (rigid\ lid)} \\ ! w &=& d\eta/dt &{\rm (old\ free\ surface)}\\ ! w &=& d\eta/dt - F_W &{\rm (variable\ thickness\ surface\ layer)} ! \end{array}\end{equation} ! ! At the surface, the time change of the surface height ! in T columns satisfies the barotropic continuity equation: ! \begin{equation} ! d\eta/dt + \nabla\cdot(H{\bf U}) = 0, ! \end{equation} ! where ${\bf U} = (U,V)$ is the barotropic velocity. ! ! On the second pass of a matsuno timestep, the predicted pressure ! from the 1st pass is contained in array PSURF(:,:,newtime) ! (PSURF(:,:,curtime) is not updated until the end of the matsuno ! step). ! ! Note that DH, DHU represents $d\eta/dt - F_W$ at ! T, U points in variable thickness surface layer ! ! !REVISION HISTORY: ! same as module ! !OUTPUT PARAMETERS: real (r8), dimension(nx_block,ny_block,max_blocks_clinic), & intent(out) :: & DH, DHU ! change in surface height (-Fw) at T,U points !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (int_kind) :: & iblock ! block index type (block) :: & this_block ! block information for current block real (r8), dimension(nx_block,ny_block) :: & WORK ! temporary work space ! real (r8), dimension(nx_block,ny_block) :: WORKX,WORKY ! real (r8) :: div_error, div_mag, dh_dt,dfw,dhtot !----------------------------------------------------------------------- ! ! calculate dh/dt or dh/dt - Fw at T points ! !----------------------------------------------------------------------- !$OMP PARALLEL DO PRIVATE(iblock, this_block, WORK) do iblock = 1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) select case (sfc_layer_type) case (sfc_layer_varthick) if (mix_pass == 2) then DH(:,:,iblock) = (PSURF(:,:,newtime,iblock) - & PSURF(:,:,curtime,iblock))/(grav*dtp) - & FW(:,:,iblock) else DH(:,:,iblock) = (PSURF(:,:,curtime,iblock) - & PSURF(:,:,oldtime,iblock))/(grav*dtp) - & FW_OLD(:,:,iblock) endif case (sfc_layer_rigid) DH(:,:,iblock) = c0 case (sfc_layer_oldfree) if (mix_pass == 2) then DH(:,:,iblock) = (PSURF(:,:,newtime,iblock) - & PSURF(:,:,curtime,iblock))/(grav*dtp) else DH(:,:,iblock) = (PSURF(:,:,curtime,iblock) - & PSURF(:,:,oldtime,iblock))/(grav*dtp) endif end select !----------------------------------------------------------------------- ! ! Some debugging diagnostics ! ...check barotropic continuity equation: ! !----------------------------------------------------------------------- ! ! WORKX = HU(:,:,iblock)*UBTROP(:,:,curtime,iblock) ! WORKY = HU(:,:,iblock)*VBTROP(:,:,curtime,iblock) ! call div(1,WORK,WORKX,WORKY,this_block) ! WORKX = (DH(:,:,iblock)*TAREA(:,:,iblock) + WORK)**2 ! WORKY = TAREA(:,:,iblock)**2 ! div_error = sqrt(global_sum(WORKX,distrb_clinic,field_loc_center,RCALCT)/& ! global_sum(WORKY,distrb_clinic,field_loc_center,RCALCT)) ! WORKX = (DH(:,:,iblock)*TAREA(:,:,iblock))**2 ! WORKY = TAREA(:,:,iblock)**2 ! div_mag=sqrt(global_sum(WORKX,distrb_clinic,field_loc_center,RCALCT)/& ! global_sum(WORKY,distrb_clinic,field_loc_center,RCALCT)) ! if (my_task == master_task) then ! write(stdout,*) ' div error: ', div_error, div_mag ! endif ! WORK = -WORK/TAREA(:,:,iblock) ! ! WORKX = TAREA(:,:,iblock)*(PSURF(:,:,curtime,iblock) - & ! PSURF(:,:,oldtime,iblock))/grav ! WORKY = dtp*TAREA(:,:,iblock)*FW_OLD(:,:,iblock) ! WORK = TAREA(:,:,iblock)*PSURF(:,:,curtime,iblock)/grav ! dh_dt = global_sum(WORKX,distrb_clinic,field_loc_center,RCALCT)/area_t ! dfw = global_sum(WORKY,distrb_clinic,field_loc_center,RCALCT)/area_t ! dhtot = global_sum(WORK ,distrb_clinic,field_loc_center,RCALCT)/area_t ! if (my_task == master_task) then ! write(stdout,*) ' ' ! write(stdout,*) 'mean surface height ', dhtot ! write(stdout,*) 'mean change in surface height ', dh_dt ! write(stdout,*) 'mean change from freshwater ', dfw ! endif ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! calculate dh/dt (or dh/dt - Fw) at U points as the area-weighted ! average of value at T points. ! !----------------------------------------------------------------------- call tgrid_to_ugrid(DHU(:,:,iblock),DH(:,:,iblock), iblock) where (.not. CALCU(:,:,iblock)) DHU(:,:,iblock) = c0 ! zero on land !----------------------------------------------------------------------- ! ! accumulate surface height tavg diagnostics if requested ! !----------------------------------------------------------------------- if (accumulate_tavg_now(tavg_SSH) ) then WORK = PSURF(:,:,curtime,iblock)/grav call accumulate_tavg_field(WORK, tavg_SSH, iblock, 1) endif if (accumulate_tavg_now(tavg_SSH2) ) then WORK = (PSURF(:,:,curtime,iblock)/grav)**2 call accumulate_tavg_field(WORK, tavg_SSH2, iblock, 1) endif if (accumulate_tavg_now(tavg_H3) ) then WORK = ((GRADPX(:,:,curtime,iblock)/grav)**2 + & (GRADPY(:,:,curtime,iblock)/grav)**2) call accumulate_tavg_field(WORK, tavg_H3, iblock, 1) endif !----------------------------------------------------------------------- ! ! accumulate surface height movie diagnostics if requested ! !----------------------------------------------------------------------- if (movie_requested(movie_SSH) .and. mix_pass /= 1) then WORK = PSURF(:,:,curtime,iblock)/grav call update_movie_field(WORK, movie_SSH, iblock, 1) endif end do ! block loop !$OMP END PARALLEL DO !----------------------------------------------------------------------- end subroutine dhdt !*********************************************************************** end module surface_hgt !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||