!-------------------------------------------------------------------------- ! Copyright (C) 2005 SPMODEL Development Group. All rights reserved. !-------------------------------------------------------------------------- ! !Title: Two-dimensional barotropic model on a double cyclic region ! ! The time integration is performed with the Adams Basthforth scheme ! (non viscous term) and Crank-Nicolson scheme (viscous term). ! The beta parameter is constant. ! Program is organized for test3 (test for modon solutions). ! This program should be complied with Fujitsu Fortran90 compiler ! and linked with SSL2 library since the subroutines for calculation ! of Bessel functions are used. ! !History: ! 2005/10/25 Shin-ichi Takehiro, created. ! ! The governing equation: ! ! d\nabla^2\psi/dt ! = -J(\psi,\nabla^2\psi) - beta d\psi/dx ! +(-1)^{p+1}\nu_{2p}\nabla^2p \nabla^2\psi ! program plbaro_beta_abcn_test3 !== Modules specification ================================================ use ee_module use gt4_history, only : GT_HISTORY, HistoryCreate, HistoryPut, HistoryGet, & HistoryAddVariable, HistoryClose, HistoryAddAttr use dc_trace, only : SetDebug, BeginSub, EndSub, DbgMessage use dc_message, only : MessageNotify implicit none !== Declaration ============================================================ !---- Variables(grid data) ---- real(8), allocatable :: yx_VelX(:,:) ! Longitudinal velocity real(8), allocatable :: yx_VelY(:,:) ! Latitudinal velocity real(8), allocatable :: yx_Vor(:,:) ! Vorticity(vertical) real(8), allocatable :: yx_StrFunc(:,:) ! Stream function real(8), allocatable :: yx_VorAnalytic(:,:) ! Analytic solution of Vorticity !---- Variables(spectral data) ---- real(8), allocatable :: ee_Vor(:,:) ! Vorticity(vertical) real(8), allocatable :: ee_StrFunc(:,:) ! Stream function real(8), allocatable :: ee_DVorDt(:,:) ! Time variation of vor.(t) real(8), allocatable :: ee_DVorDtB(:,:) ! Time variation of vor.(t-delta_t) real(8), allocatable :: ee_VorAnalytic(:,:) ! Analytic solution of Vorticity real(8), allocatable :: ee_VorInit(:,:) ! Initial value of Vorticity !---- Fixed parameters ----- real(8), parameter :: pi = 3.141592653589793D0 ! Circle ratio character(len=20) :: DbgMessageFmt='*** DbgMESSAGE ***' real(8), parameter :: vmiss = -999.0 ! Missing value !---- NAMELIST variables ---- ! Filename for NAMELIST reading (should not be changed) character(len=30),parameter :: nmlfile='plbaro-beta_abcn_test3.nml' logical :: Verbose=.false. ! Verbose message On/Off logical :: DebugOn=.false. ! Debug message On/Off namelist /message/ Verbose, DebugOn ! !-- Grid, Spectral -- integer :: km=21 ! Truncation wavenumber(X) integer :: lm=21 ! Truncation wavenumber(Y) integer :: im=64 ! Number of grid points (X) (>3*km) integer :: jm=64 ! Number of grid points (>3*lm) namelist /gridset/ km, lm, im, jm !-- Physical parameters -- real(8) :: XLength=1.0D0 ! Domain length (X) real(8) :: YLength=1.0D0 ! Domain length (Y) real(8) :: Beta=0.0D0 ! Beta parameter integer :: HVOrder=1 ! Order of hyper viscosity ! (1 means normal viscosity, ! Order of hor. Laplacian) real(8) :: HVisc=1.0D0 ! Hyper visc. coefficient namelist /physics/ XLength, YLength, Beta, HVOrder, HVisc ! -- Initial condition -- character(len=100) :: initial_file='' ! Filename of initial data ! (They are calculated if empty) real :: initial_time=0.0 ! Initial time namelist /initial/ initial_file, initial_time ! ! -- Time integration -- real(8) :: delta_t=1.0e-7 ! Time step interval integer :: nstep=2000 ! Total time step namelist /tint/ delta_t, nstep ! ! -- History output -- character(len=100) :: hst_file= '' ! Filename of history output character(len=100) :: title = & ! Title 'Test of linear terms of 2-dim barotropic model on a double-cyclic domain' integer :: hst_intstep=200 ! Output interval namelist /history/ hst_file, title, hst_intstep character(len=100) :: rst_file='' ! Filename for restart integer :: rst_intstep=200 ! Output interval namelist /restart/ rst_file, rst_intstep !--Initial data: one modon real(8) :: Radius = 1.0 ! Radius of Modon real(8) :: Speed = 1.0 ! propagation speed of Modon real(8) :: XCenter = 0.5 ! Central position(X) real(8) :: YCenter = 0.5 ! Central position(Y) namelist /modon/ Radius, Speed, XCenter, YCenter !---- Other variables ---- real(8), allocatable :: ee_HVisc(:,:) ! Hyper visc. coefficient integer :: it=0 ! Time step real(8) :: time ! Time in the model integer :: k, l ! Wavenumbers integer :: i, j ! Index of grid points type(GT_HISTORY) :: hst_rst ! Restart GT_HISTORY variable !---------------- Reading NAMELIST --------------------- write(6,nml=message) open(10,file=nmlfile,status='OLD') read(10,nml=message) ; write(6,nml=message) ; close(10) if (verbose) write(6,nml=gridset) open(10,file=nmlfile,status='OLD') read(10,nml=gridset) ; write(6,nml=gridset) ; close(10) if (verbose) write(6,nml=physics) open(10,file=nmlfile,status='OLD') read(10,nml=physics) ; write(6,nml=physics) ; close(10) if (verbose) write(6,nml=initial) open(10,file=nmlfile,status='OLD') read(10,nml=initial) ; write(6,nml=initial) ; close(10) if (verbose) write(6,nml=tint) open(10,file=nmlfile,status='OLD') read(10,nml=tint) ; write(6,nml=tint) ; close(10) if (verbose) write(6,nml=history) open(10,file=nmlfile,status='OLD') read(10,nml=history) ; write(6,nml=history) ; close(10) if (verbose) write(6,nml=restart) open(10,file=nmlfile,status='OLD') read(10,nml=restart) ; write(6,nml=restart) ; close(10) if (verbose) write(6,nml=modon) open(10,file=nmlfile,status='OLD') read(10,nml=modon) ; write(6,nml=modon) ; close(10) !---------------- Debugging option ----------------- if (DebugOn) then call SetDebug end if !------------------ Allocate variables --------------------- allocate(yx_VelX(0:jm-1,0:im-1),yx_VelY(0:jm-1,0:im-1)) allocate(yx_Vor(0:jm-1,0:im-1),yx_StrFunc(0:jm-1,0:im-1)) allocate(yx_VorAnalytic(0:jm-1,0:im-1)) allocate(ee_Vor(-lm:lm,-km:km),ee_StrFunc(-lm:lm,-km:km)) allocate(ee_DVorDt(-lm:lm,-km:km),ee_DVorDtB(-lm:lm,-km:km)) allocate(ee_HVisc(-lm:lm,-km:km)) allocate(ee_VorAnalytic(-lm:lm,-km:km),ee_VorInit(-lm:lm,-km:km)) !------------------ Coordinates ----------------------- call DbgMessage(fmt='call %c', c1='ee_initial') call ee_Initial(im,jm,km,lm,0.0D0,XLength,0.0D0,YLength) !------------------ Physical parameters ----------------------- ee_Vor = 1.0D0 ee_HVisc = (-1)**HVOrder * HVisc*(ee_Lapla_ee(ee_Vor))**HVOrder !------------------- Initial condition ---------------------- time = initial_time if ( initial_file == "") then ! ee_Vor, ee_DVorDtB are given internally if filename is empty call set_initial_values else ! Set initial value from the restart file call HistoryGet( trim(initial_file), 'ee_vor', ee_Vor, time ) call HistoryGet( trim(initial_file), 'ee_dvordt_before',ee_DVorDtB,time ) endif ee_Strfunc = ee_LaplaInv_ee(ee_Vor) !----------- Time integration(Leap frog + Crank Nicolson scheme) ---------- call output_restart_init call output_history_init if ( initial_file == '' ) call output_history ! Output if initial values ! are given internally call DbgMessage(fmt='%c %c', & & c1=DbgMessageFmt, & & c2='Time integration starts.') do it=1,nstep time = initial_time + it * delta_t ! Evaluate time variation of vorticity ee_DVorDt = - ee_Jacobian_ee_ee(ee_StrFunc,ee_Vor) & - Beta * ee_Dx_ee(ee_StrFunc) ee_Vor = ( (1.0D0 - ee_HVisc*delta_t/2)*ee_Vor & + 3.0D0 * delta_t/2.0D0 * ee_DVorDt & - 1.0D0 * delta_t/2.0D0 * ee_DVorDtB ) & /(1.0D0 + ee_HVisc*delta_t/2) ! Transposition of the variable ee_DVorDtB = ee_DVorDt ! Calculating stream function from vorticity ee_StrFunc = ee_LaplaInv_ee(ee_Vor) if(mod(it,hst_intstep) .eq. 0)then ! History output call output_history endif if(mod(it,rst_intstep) .eq. 0)then ! Restart output call output_restart endif enddo call DbgMessage(fmt='%c %c', & & c1=DbgMessageFmt, & & c2='Time integration end.') if(.not. mod(it-1,rst_intstep) .eq. 0)then ! Final output call output_restart endif call output_restart_close call output_history_close ! main program end !----------------------------------------------------------------------------- ! subprograms below contains !=========================== Modon solution ============================ ! ! Vorticity of modon solution ! subroutine set_VorMod( & yx_VorMod, & !(out) Vorticity of Modon Ra, & !(in) Radius of Modon Cs, & !(in) propagation speed of Modon Xc, & !(in) Central position(X) Yc & !(in) Central position(Y) ) real(8), intent(out):: yx_VorMod(0:,0:) ! Vorticity of Modon real(8), intent(in) :: Ra ! Radius of Modon real(8), intent(in) :: Cs ! propagation speed of Modon real(8), intent(in) :: Xc ! Central position(X) real(8), intent(in) :: Yc ! Central position(Y) ! Internal parameter real(8),parameter :: eps=1.0D-10 ! Termination error for ! determination of K ! Working variables real(8) :: Kp ! A parameter of the modon solution real(8) :: Qp ! A parameter of the modon solution real(8) :: J1KR, J1K ! 1st kind Bessel functiono real(8) :: K1KR, K1K ! 2nd kind modified Bessel function real(8) :: Rad ! Radius from the center real(8) :: SinTheta ! Plane polar coordinate, \sin\theta integer :: ICON ! Error code of subroutines real(8) :: DX, DY ! Working areas Qp = Ra*sqrt(Beta/Cs) Kp = Modon_Kparameter(Qp,eps) call dbj1(Kp,J1K,ICON) call dbk1(Qp,K1K,ICON) do i=0,im-1 do j=0,jm-1 if ( yx_X(j,i)-mod(Xc,XLength) .LT. -XLength/2 )then DX=yx_X(j,i)+XLength-mod(Xc,XLength) else if ( yx_X(j,i)-mod(Xc,XLength) .GT. XLength/2 )then DX=yx_X(j,i)-XLength-mod(Xc,XLength) else DX=yx_X(j,i)-mod(Xc,XLength) endif if ( yx_Y(j,i)-mod(Yc,YLength) .LT. -YLength/2 )then DY=yx_Y(j,i)+YLength-mod(Yc,YLength) else if ( yx_Y(j,i)-mod(Yc,YLength) .GT. YLength/2 )then DY=yx_Y(j,i)-YLength-mod(Yc,YLength) else DY=yx_Y(j,i)-mod(Yc,YLength) endif Rad = sqrt( DX**2 + DY**2 ) if ( Rad .NE. 0.0D0 ) then SinTheta = DY/Rad else SinTheta = 0.0D0 endif if ( Rad .LE. Ra ) then call dbj1(Kp*Rad/Ra,J1KR,ICON) yx_VorMod(j,i) = - Cs * Qp**2/Ra * SinTheta * J1KR/J1K else call dbk1(Qp*Rad/Ra,K1KR,ICON) yx_VorMod(j,i) = - Cs * Qp**2/Ra * SinTheta * K1KR/K1K endif enddo enddo end subroutine set_VorMod ! ! Solving Modon dispersion relation. S.Takehiro Oct, 21, 2005. ! ! Calculating parameter k by giving the value of parameter q ! with the relation, ! -J_2(k)/kJ_1(k) = K_2(q)/qK_1(q). ! The bi-section method is used in the function. ! function Modon_Kparameter(Qp,eps) real(8), intent(in) :: Qp ! Parameter q real(8), intent(in) :: eps ! Termination error real(8) :: Modon_Kparameter ! Parameter k ! Working variables real(8) :: k1 = 3.83171 ! Initial value: J_1 1st zero point real(8) :: k2 = 5.13562 ! Initial value: J_2 1st zero point real(8) :: kp ! Guess of the parameter real(8) :: J1K, J2K ! 1st kind Bessel function real(8) :: K1Q, K2Q ! 2nd kind modified Bessel function integer :: ICON ! Error code of subroutines real(8) :: f1, f2, fp ! Working areas 1000 continue call dbk1(Qp,K1Q,ICON) call dbkn(Qp,2,K2Q,ICON) call dbj1(k1,J1K,ICON) call dbjn(k1,2,J2K,ICON) f1 = J2K + K2Q/(Qp*K1Q)*k1*J1K ! J_2(k) + K_2(q)/qK_1(q) * kJ_1(k) =0 call dbj1(k2,J1K,ICON) call dbjn(k2,2,J2K,ICON) f2 = J2K + K2Q/(Qp*K1Q)*k2*J1K ! J_2(k) + K_2(q)/qK_1(q) * kJ_1(k) =0 if ( f1*f2 .gt. 0.0D0 ) then write(6,*) 'f1 and f2 are the same sign.' stop endif ! kp = (abs(f1)*k2 + abs(f2)*k1)/abs(f1-f2) kp = (k1+k2)/2 call dbj1(kp,J1K,ICON) call dbjn(kp,2,J2K,ICON) fp = J2K + K2Q/(Qp*K1Q)*kp*J1K ! J_2(k) + K_2(q)/qK_1(q) * kJ_1(k) =0 if ( abs(k1-k2)/abs(k1) .lt. eps ) then Modon_Kparameter=kp return end if if ( f1*fp .gt. 0.0D0 ) then k1 = kp else k2 = kp endif goto 1000 end function Modon_Kparameter !=========================== Initial condition ============================ ! ! Default initial values when the restart file is not specified ! subroutine set_initial_values ! gives w_Vor call set_VorMod( & yx_Vor, & !(out) Vorticity of Modon Radius, & !(in) Radius of Modon Speed, & !(in) Propagation speed of Modon XCenter, & !(in) Central position(X) YCenter & !(in) Central position(Y) ) ee_Vor = ee_yx(yx_Vor) ee_DVorDtB = - ee_Jacobian_ee_ee(ee_StrFunc,ee_Vor) & - Beta * ee_Dx_ee(ee_StrFunc) end subroutine set_initial_values !=========================== Output for restart ============================ ! ! Initialization of restart output ! subroutine output_restart_init call HistoryCreate( & file=trim(rst_file), & title=trim(title), & source='plbaro-beta_abcn_test3.f90 (2005/10/25)', & institution='GFD_Dennou Club SPMODEL project', & dims=(/'x','y','k','l','t'/), & dimsizes=(/im,jm,2*km+1,2*lm+1,0/),& longnames=(/'X ','Y ',& 'X-wavenumber','Y-wavenumber',& 'time '/),& units=(/'1','1','1','1','1'/), & origin=real(time), interval=real(rst_intstep*delta_t), & xtypes=(/'real'/), history=hst_rst) !---- Define, output coordinate variables ---- call HistoryPut('x',x_X, hst_rst) ! Output variable call HistoryAddattr('x','topology','circular', hst_rst) ! Cyclic coordinate call HistoryAddattr('x','modulo',XLength, hst_rst) ! Cyclic coordinate call HistoryPut('y',y_Y, hst_rst) ! Output variable call HistoryAddattr('y','topology','circular', hst_rst) ! Cyclic coordinate call HistoryAddattr('y','modulo',YLength, hst_rst) ! Cyclic coordinate call HistoryPut('k',(/(dble(k),k=-2*km,2*km)/), hst_rst)! Output variable call HistoryPut('k',(/(dble(l),l=-2*lm,2*lm)/), hst_rst)! Output variable call HistoryAddVariable( & ! Define variable varname='x_weight', dims=(/'x'/), & longname='Weight for integration in X', & units='1', xtype='double',history=hst_rst) call HistoryPut('x_weight',x_X_Weight,hst_rst) ! Output variable call HistoryAddVariable( & ! Define variable varname='y_weight', dims=(/'y'/), & longname='Weight for integration in Y', & units='1', xtype='double',history=hst_rst) call HistoryPut('y_weight',y_Y_Weight,hst_rst) ! Output variable !---- Define physical variables ---- call HistoryAddVariable( & ! Define variable varname='ee_vor', dims=(/'l','k','t'/), & longname='Vorticity', & units='1', xtype='double', history=hst_rst) call HistoryAddVariable( & ! Define variable varname='ee_dvordt_before', dims=(/'l','k','t'/), & longname='Time variation of Vorticity (1step before)', & units='1', xtype='double', history=hst_rst) !---- Output of exp. parameters as global attribtes ---- call HistoryAddAttr('x','+delta_t', delta_t ,hst_rst) call HistoryAddAttr('x','+Beta', Beta ,hst_rst) call HistoryAddAttr('x','+HVOrder', HVOrder ,hst_rst) call HistoryAddAttr('x','+HVisc', HVisc ,hst_rst) end subroutine output_restart_init ! ! Restart output ! subroutine output_restart write(6,*) ' Restart file output at it = ',it, ' time = ', time call HistoryPut('t',real(time),hst_rst) !---- Output physical variables ---- call HistoryPut('ee_vor', ee_Vor, hst_rst) call HistoryPut('ee_dvordt_before', ee_DVorDtB, hst_rst) end subroutine output_restart ! ! Closing restart output ! subroutine output_restart_close call HistoryClose(hst_rst) end subroutine output_restart_close !=========================== History output ============================ ! ! Initialization of history output ! subroutine output_history_init call HistoryCreate( & file=trim(hst_file), & title=trim(title), & source='plbaro-beta_abcn_test3.f90 (2005/10/25)', & institution='GFD_Dennou Club SPMODEL project', & dims=(/'x','y','k','l','t'/), & dimsizes=(/im,jm,2*km+1,2*lm+1,0/),& longnames=(/'X ','Y ',& 'X-wavenumber','Y-wavenumber',& 'time '/),& units=(/'1','1','1','1','1'/), & origin=real(time), interval=real(rst_intstep*delta_t), & xtypes=(/'real'/)) !---- Define, output coordinate variables ---- call HistoryPut('x',x_X ) ! Output variable call HistoryAddattr('x','topology','circular') ! Cyclic coordinate call HistoryAddattr('x','modulo',XLength ) ! Cyclic coordinate call HistoryPut('y',y_Y ) ! Output variable call HistoryAddattr('y','topology','circular') ! Cyclic coordinate call HistoryAddattr('y','modulo',YLength ) ! Cyclic coordinate call HistoryPut('k',(/(dble(k),k=-2*km,2*km)/) ) ! Output variable call HistoryPut('k',(/(dble(l),l=-2*lm,2*lm)/) ) ! Output variable call HistoryAddVariable( & ! Define variable varname='x_weight', dims=(/'x'/), & longname='Weight for integration in X', & units='1', xtype='double') call HistoryPut('x_weight',x_X_Weight) ! Output variable call HistoryAddVariable( & ! Define variable varname='y_weight', dims=(/'y'/), & longname='Weight for integration in Y', & units='1', xtype='double') call HistoryPut('y_weight',y_Y_Weight) ! Output variable !---- Define, output parameters ---- call HistoryAddVariable( & ! Define variable varname='hvisc', dims=(/'l','k'/), & longname='hyper diffusivity', units='1', xtype='double') call HistoryPut('hvisc',ee_HVisc) ! Output variable !---- Define physical variables ---- call HistoryAddVariable( & ! Define variable varname='vor', dims=(/'x','y','t'/), & longname='Vorticity', units='1', xtype='double') call HistoryAddVariable( & ! Define variable varname='strfunc', dims=(/'x','y','t'/), & longname='Stream function', units='1', xtype='double') call HistoryAddVariable( & ! Define variable varname='velx', dims=(/'x','y','t'/), & longname='x-velocity', units='1', xtype='double') call HistoryAddVariable( & ! Define variable varname='vely', dims=(/'x','y','t'/), & longname='y-velocity', units='1', xtype='double') call HistoryAddVariable( & ! Define variable varname='vor_analytic', dims=(/'x','y','t'/), & longname='Vorticity(analytic)', units='1', xtype='double') !---- Define diagnostic variables ---- call HistoryAddVariable( & ! Define variable varname='ek', dims=(/'t'/), & longname='mean kinetic energy', units='1', xtype='double') call HistoryAddVariable( & ! Define variable varname='ens', dims=(/'t'/), & longname='mean enstrophy', units='1', xtype='double') call HistoryAddVariable( & ! Define variable varname='ektot', dims=(/'t'/), & longname='total kinetic energy', units='1', xtype='double') call HistoryAddVariable( & ! Define variable varname='enstot', dims=(/'t'/), & longname='total enstrophy', units='1', xtype='double') call HistoryAddVariable( & ! Define variable varname='l2vor', dims=(/'t'/), & longname='vorrticity error (type2)', & units='1', xtype='double') !---- Output of exp. parameters as global attribtes ---- call HistoryAddAttr('x','+delta_t', delta_t ) call HistoryAddAttr('x','+Beta', Beta ) call HistoryAddAttr('x','+HVOrder', HVOrder ) call HistoryAddAttr('x','+HVisc', HVisc ) end subroutine output_history_init ! ! History output ! subroutine output_history write(6,*) ' History file output at it = ',it, ' time = ', time call HistoryPut('t',real(time)) !---- Output physical variables ---- call set_VorMod( & yx_VorAnalytic, & !(out) Vorticity of Modon Radius, & !(in) Radius of Modon Speed, & !(in) Propagation speed of Modon XCenter+Speed*time, & !(in) Central position(X) YCenter & !(in) Central position(Y) ) yx_VorAnalytic = yx_ee(ee_yx(yx_VorAnalytic)) yx_Vor = yx_ee(ee_Vor) yx_VelX = -yx_ee(ee_Dy_ee(ee_StrFunc)) yx_VelY = yx_ee(ee_Dx_ee(ee_StrFunc)) yx_StrFunc = yx_ee(ee_StrFunc) call HistoryPut('velx',transpose(yx_VelX)) call HistoryPut('vely',transpose(yx_VelY)) call HistoryPut('vor',transpose(yx_Vor)) call HistoryPut('strfunc',transpose(yx_StrFunc)) call HistoryPut('vor_analytic',transpose(yx_VorAnalytic)) !---- Output diagnostic variables ---- call HistoryPut('ek', AvrYX_yx((yx_VelX**2+yx_VelY**2)/2.0d0)) call HistoryPut('ens', AvrYX_yx(yx_Vor**2/2.0d0)) call HistoryPut('ektot', IntYX_yx((yx_VelX**2+yx_VelY**2)/2.0d0)) call HistoryPut('enstot', IntYX_yx(yx_Vor**2/2.0d0)) call HistoryPut('l2vor', & sqrt(AvrYX_yx((yx_Vor-yx_VorAnalytic)**2)& /AvrYX_yx(yx_VorAnalytic**2))) end subroutine output_history ! ! Closing history output ! subroutine output_history_close call HistoryClose end subroutine output_history_close end program plbaro_beta_abcn_test3