!---------------------------------------------------------------------- ! Copyright (c) 2002--2005 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- ! Sample program for SPML(based on gt4f90io and ISPACK) ! ! 2002/08/19 S.Takehiro ! 2004/01/26 M.Odaka ! 2004/06/17 M.Odaka ! 2005/03/16 S.Takehiro ! ! Solving 2-D Boussinesq fluid system (FF problem) ! d\zeta/dt + J(\psi,\zeta) = PrRa dT/dx + Pr\nabla\zeta ! dT/dt + J(\psi,T) - d\psi/dx = \nabla T ! \nabla\psi = \zeta ! psi = zeta = dT/dy = 0 at y=0,1 ! program bsncnv2 use esc_module use gt4_history implicit none !---- Spartial resolution ---- integer, parameter :: km=32, lm=10 ! Truncation wavenumbers(X,Y) integer, parameter :: im=128, jm=16 ! Grid points(X,Y) !---- Variables ---- real(8) :: yx_Psi(0:jm,0:im-1) ! Grid (Stream function) real(8) :: yx_Temp(0:jm,0:im-1) ! Grid (Temperature disturbance) real(8) :: yx_Zeta(0:jm,0:im-1) ! Grid (Vorticity) real(8) :: es_PsiA(-km:km,lm) ! Spectral (Stream function,t) real(8) :: ec_TempA(-km:km,0:lm)! Spectral (Temp. disturb.,t) real(8) :: es_ZetaA(-km:km,lm) ! Spectral (Vorticity,t) real(8) :: es_PsiB(-km:km,lm) ! Spectral (Stream function,t+1) real(8) :: ec_TempB(-km:km,0:lm)! Spectral (Temp. disturb.,t+1) real(8) :: es_ZetaB(-km:km,lm) ! Spectral (Vorticity,t+1) !---- Parameters for coodinates ---- real(8), parameter :: xmin=0.0, xmax=8.0 ! Domain length (X) real(8), parameter :: ymin=0.0, ymax=1.0 ! Domain length (Y) !---- Parameters for time integration ---- real(8), parameter :: dt=1e-3 ! Time step interval integer, parameter :: nt=20000, ndisp=500 ! Total, display time step !---- Physical parameters ---- real(8), parameter :: Ra=1.0e4 ! Rayleigh number real(8), parameter :: Pr=1.0 ! Prandtl number integer :: it ! DO variable !---------------- Setting coodinates --------------------- call esc_initial(im,jm,km,lm,xmin,xmax,ymin,ymax) ! Initialize module !-------------- Setting initial values ------------------ yx_Temp = 0.0 ; yx_Temp(jm/2,im/2) = 0.01 ! temperature perturbation yx_Psi = 0.0 yx_Zeta = 0.0 ec_TempA = ec_yx(yx_Temp) ; ec_TempB = ec_TempA es_PsiA = es_yx(yx_Psi) ; es_PsiB = es_PsiA es_ZetaA = es_yx(yx_Zeta) ; es_ZetaB = es_ZetaA call output_gtool4_init ! Initialize history call output_gtool4 ! Output initial data !------------------- Time integration ---------------------- do it=1,nt ec_TempA = ec_TempB + & dt*( -ec_Jacobian_es_ec(es_PsiB,ec_TempB) & + ec_yx(yx_es(es_dx_es(es_PsiB))) & + ec_Lapla_ec(ec_TempB) ) es_ZetaA = es_ZetaB + & dt*( - es_Jacobian_es_es(es_PsiB,es_ZetaB) & + Pr * Ra * es_yx(yx_ec(ec_Dx_ec(ec_TempB))) & + Pr * es_Lapla_es(es_ZetaB) ) es_PsiA = es_LaplaInv_es(es_ZetaA) ec_TempB = ec_TempA es_PsiB = es_PsiA es_ZetaB = es_ZetaA if(mod(it,ndisp) .eq. 0)then ! Output call output_gtool4 endif enddo call output_gtool4_close stop contains subroutine output_gtool4_init call HistoryCreate( & ! History create file='bsncnv-ff-1.nc', title='Boussinesq convection (FF)', & source='Sample program of gtool_history/gtool4', & institution='GFD_Dennou Club davis/spmodel project',& dims=(/'x','y','t'/), dimsizes=(/im,jm+1,0/), & longnames=(/'X-coordinate','Y-coordinate','time '/),& units=(/'1','1','1'/), & origin=0.0, interval=real(ndisp*dt) ) call HistoryPut('x',x_X) ! Output variable call HistoryAddattr('x','topology','circular') ! Cyclic coordinate call HistoryAddattr('x','modulo',xmax-xmin) ! Cyclic coordinate call HistoryPut('y',y_Y) ! Output variable call HistoryAddVariable( & ! Define variable varname='psi', dims=(/'x','y','t'/), & longname='stream function', units='1', xtype='double') call HistoryAddVariable( & ! Define variable varname='zeta', dims=(/'x','y','t'/), & longname='vorticity', units='1', xtype='double') call HistoryAddVariable( & ! Define variable varname='temp', dims=(/'x','y','t'/), & longname='temperature', units='1', xtype='double') end subroutine output_gtool4_init subroutine output_gtool4 yx_Temp = yx_ec(ec_TempA) + 1 - yx_Y ! Output is total temperature yx_Psi = yx_es(es_PsiA) yx_Zeta = yx_es(es_ZetaA) write(6,*) 'it = ',it call HistoryPut('t',real(it*dt)) call HistoryPut('psi',transpose(yx_Psi)) call HistoryPut('zeta',transpose(yx_Zeta)) call HistoryPut('temp',transpose(yx_Temp)) end subroutine output_gtool4 subroutine output_gtool4_close call HistoryClose end subroutine output_gtool4_close end program bsncnv2