C This file serves for development of numerical solution of parametric C C description of mean curvature flow C C Prepis: C 1. zmena rozmeru poli C 2. uprava funkce prave strany C#######################################################################C C C C Method of lines for solving the graph variant of the elastic C C flow in the form: C C C C alpha v_Gamma_Phi = Delta_Gamma_Phi ( Kappa_Gamma_Phi + F ) C C C C Version 0.2.2.beta C C Playing with the i.c. shape C C-----------------------------------------------------------------------C C Comes from the isotropic versions of: C C Novebre 6, 2001, April 18, 2002, April 26, 2002, October 28, 2002 C C May 29, 2003, August 26, 2005, September 19-22, 2005 C C Date: June 6, 2006, September 4, 2006 C C C C#######################################################################C C C C Dr.Ing. Michal Benes C C Dept of Maths, FJFI CVUT C C (C) 2001 - 2006 C C C C#######################################################################C C f90 -K -O3 +U77 -o elast22_f_hp elast22_f.f C C ifc -pc80 -FI -o elast22_f_i elast22_f.f C C ifc -pc80 -FI -w95 -Vaxlib -o elast22_f_i elast22_f.f C C ifc -prec_div -FI -w95 -Vaxlib -o elast22_f_i elast22_f.f C C ifc -mp -FI -w95 -tpp7 -save -Vaxlib -o elast22_f_i elast22_f.f C C f90 -K -O3 +Oall +Ovectorize +DD64 +U77 -o C C elast22_f_ani_hp elast22_f_ani.f C PROGRAM Anisotropic_Elastic_by_Lines C time measurement CHARACTER*8 time_string INTEGER*2 hour, minute, second, hundrd C date and time CHARACTER*24 FDATE REAL*4 ETIME, ELAPSED, SYS(2) INTEGER date(8) CHARACTER (LEN=10) BEN(3) C simulation settings DOUBLE PRECISION T_STEP(500000), check_sum INTEGER*4 N_TST, TST_N DOUBLE PRECISION H1, H2, DT, F, ZETA, THETA_0, R, EPS, L1, L2, / Phi(-1:401,-1:401), DELTA, STEP, ALPHA, BLOW, / SCALE, a, b, c C pole Phi bude mit rozmer 2 x -1:401??? INTEGER FRQ, FIN_ID, WRITE_ID, CUR_ID, FLD_ID, CHK_ID, / SAVE_ID, LOAD_ID, SYS_ID LOGICAL FIRST COMMON /NUMER/ H1, H2, DT, N1, N2, EPS /MERSN/ STEP /CHECK/ BLOW / /PHYS/ ALPHA, F, ZETA, THETA_0, M_FOLD, R / /PATTERN/ SCALE /ID/ WRITE_ID /STORE/ SAVE_ID, LOAD_ID / /SYS_DATA/ SYS_ID, N_TST, T_STEP, TST_N, check_sum / /ANIS/ a, b, c C-----------------------------------------------------------------------C WRITE_ID = 1 WRITE (*,*) 'Method-of-lines scheme for solving' WRITE (*,*) 'the graph formulation of the elastic flow' WRITE (*,*) 'Version 0.2.2.beta' WRITE (*,*) 'Reading data' OPEN (96,FILE='elast22.dta',STATUS='OLD') READ (96,'(////////////2x,1F12.5,I12,5F12.5)') / F,M_FOLD,ZETA,THETA_0,R,SCALE READ (96,'(///2x,6F12.5)') L1, L2, ALPHA, a, b, c READ (96,'(//////2x,2I12,F12.5,I12,F12.5)') N1,N2,DT,NT,EPS READ (96,'(///2x,3F12.5)') DELTA, STEP, BLOW READ (96,'(//////2x,6I12///2X,3I12)') FRQ, FIN_ID, INI_ID, CUR_ID, / FLD_ID, CHK_ID, SAVE_ID, LOAD_ID, SYS_ID CLOSE(96) WRITE (*,*) 'Data loaded.' IF ((N1.GT.800).OR.(N2.GT.800)) THEN WRITE (*,*) 'Grid size exceeded!' STOP ENDIF H1 = L1/N1 H2 = L2/N2 C krok site je H1 a H2, u krivek se uzije jen jedno z nich PI = 3.14159D00 if ((STEP.gt.DT).or.(STEP.le.0.0D00)) STEP = DT C-----------------------------------------------------------------------C C timing block C system data initialization C N_TST = 0 TST_N = 0 check_sum = 0.D00 OPEN (96,FILE='elast22.dta',STATUS='OLD',ACCESS='APPEND') C PC call GetTim (hour, minute, second, hundrd) C Lin CALL DATE_AND_TIME (BEN(1),BEN(2),BEN(3),DATE) WRITE (96,*) 'Process started on:' C Lin ,DATE(3),'.',DATE(2),'.',DATE(1) ELAPSED = ETIME(SYS) WRITE (96,*)'Initial time: ' C Lin,DATE(5),':',DATE(6),':',DATE(7),'.', C Lin / DATE(8) WRITE (96,*) FDATE() write (96,*) 'Elapsed time = ', ELAPSED, ' User time = ', / SYS(1), ' System time = ', SYS(2) CLOSE (96) C initial condition C bud definuje p.p. funkci Initial_Level, nebo nahraje z dat (STORAGE) WRITE (*,*) 'Initialization in progress...' IF (LOAD_ID.NE.1) THEN CALL INITIAL_LEVEL (Phi) ELSE CALL STORAGE (Phi) ENDIF K = 0 IF (INI_ID.EQ.1) THEN C OPEN (77,FILE='field.ini',STATUS='unknown') CALL SAVE_FIELD (Phi, K) C CLOSE (77) ENDIF WRITE (*,*) 'Initialization terminated.' C----------------------------------------------------------------------C C time loop FIRST = .TRUE. DO K=1,NT CALL UPDATE_TIME_LEVEL((K-1)*DT,Phi,K*DT,N1,N2,DELTA,FIRST) IF (MOD(K,FRQ).EQ.0) THEN WRITE (*,*) 'Time level No.: ', K C ulozi Phi do souboru IF (FLD_ID.EQ.1) then CALL SAVE_FIELD (Phi, K) ENDIF ENDIF IF (SYS_ID.EQ.1) THEN OPEN(65,FILE='elast22.tst',STATUS='UNKNOWN',ACCESS='APPEND') WRITE (65,'(19H# Time level number,I10)') K DO I=1,N_TST check_sum = check_sum + T_STEP(I) WRITE(65,'(1X,I6,2X,F20.10)') TST_N+I, T_STEP(I) ENDDO CLOSE (65) WRITE (*,*) 'Partial number of time steps: ', N_TST TST_N = TST_N + N_TST N_TST = 0 ENDIF ENDDO C-----------------------------------------------------------------------C IF (FIN_ID.EQ.1) THEN WRITE (*,*) 'Saving the final form of the field Phi' C OPEN (77,FILE='field.fin',STATUS='unknown') CALL SAVE_FIELD (Phi, NT) C CLOSE (77) ENDIF IF (SAVE_ID.EQ.1) THEN CALL STORAGE (Phi) ENDIF C-----------------------------------------------------------------------C OPEN (96,FILE='elast22.dta',STATUS='OLD',ACCESS='APPEND') C PC call GetTim (hour, minute, second, hundrd) C Lin CALL DATE_AND_TIME (BEN(1),BEN(2),BEN(3),DATE) C PC WRITE (96,*)'Final time: ',hour,':',minute,':',second,'.',hundrd WRITE (96,*) 'Process terminated on:' C Lin ,DATE(3),'.',DATE(2),'.', C Lin / DATE(1) ELAPSED = ETIME(SYS) WRITE (96,*) FDATE() WRITE (96,*) 'Elapsed time = ', ELAPSED, ' User time = ', / SYS(1), ' System time = ', SYS(2) C Lin WRITE (96,*)'Final time: ' C Lin,DATE(5),':',DATE(6),':',DATE(7),'.', C Lin / DATE(8) if (SYS_ID.EQ.1) then WRITE (96,*) 'Number of time steps: ', TST_N WRITE (96,*) 'Number spatial nodes: ', / (N1+1)*(N2+1) WRITE (96,'(37H Total number of degrees of freedom: ,F25.0)') / DBLE(N1+1)*DBLE(N2+1)*DBLE(TST_N) Write (96,*) 'Control sum of time steps: ', check_sum endif CLOSE (96) WRITE (*,*) 'Program is successfully going to terminate.' WRITE (*,*) 'Time-step size at the end of computation: ', STEP STOP END C-----------------------------------------------------------------------C C zde staci prepsat pouze rozmery poli SUBROUTINE UPDATE_TIME_LEVEL (TB, Y, TE, N1, N2, DELTA, FIRST) C Runge-Kutta fourth-order scheme with automatic time stepping C DOUBLE PRECISION Y(-1:401,-1:401),Y0(0:400,0:400),K1(0:400,0:400), /K3(0:400,0:400),K4(0:400,0:400),K5(0:400,0:400),EPS,DELTA,Q DOUBLE PRECISION STEP,T,TB,TE,T0,H,H03, T_STEP(500000) INTEGER SYS_ID, N_TST LOGICAL LAST, FIRST COMMON /MERSN/ STEP COMMON /STACK/ Y0,K1,K3,K4,K5 COMMON /SYS_DATA/ SYS_ID, N_TST, T_STEP T = TB IF (FIRST) GO TO 101 H = STEP 102 IF (DABS(T-TE).LE.DABS(H)) GO TO 101 LAST = .FALSE. GO TO 103 101 LAST = .TRUE. H = TE - T 103 T0 = T DO 104 I=0,N1 DO 104 J=0,N2 104 Y0(I,J) = Y(I,J) H03 = H/3.D00 C write(*,*) 'Writing position: T=', T, ' H=', H CALL RIGHT_HAND_SIDE(T, Y, K1) DO 105 I=0,N1 DO 105 J=0,N2 105 Y(I,J) = K1(I,J)*H03 + Y0(I,J) T = T0 + H03 CALL RIGHT_HAND_SIDE(T, Y, K3) DO 106 I=0,N1 DO 106 J=0,N2 106 Y(I,J) = (K1(I,J) + K3(I,J))/2.D00*H03 + Y0(I,J) CALL RIGHT_HAND_SIDE(T, Y, K3) DO 107 I=0,N1 DO 107 J=0,N2 107 Y(I,J) = (K1(I,J)*0.375D00+K3(I,J)*1.125D00)*H03+Y0(I,J) T = T0 + H/2.D00 CALL RIGHT_HAND_SIDE(T, Y, K4) DO 108 I=0,N1 DO 108 J=0,N2 108 Y(I,J) = (K1(I,J)*1.5D00-K3(I,J)*4.5D00+K4(I,J)*6.D00)*H03+Y0(I,J) T = T0 + H CALL RIGHT_HAND_SIDE(T, Y, K5) DO 109 I=0,N1 DO 109 J=0,N2 109 Y(I,J) = ((K1(I,J)+K5(I,J))/2.D00+K4(I,J)*2.D00)*H03+Y0(I,J) EPS = 0.D00 C WRITE (*,*) 'EPS=0 ', EPS DO 110 I=0,N1 DO 110 J=0,N2 Q = DABS(K1(I,J)*0.2D00-K3(I,J)*0.9D00+K4(I,J)*0.8D00- / K5(I,J)*0.1D00) 110 IF (Q.GT.EPS) EPS = Q C write (*,*) 'Eps = ', EPS, ' Q= ', Q C write (*,*) 'Writing EPS= ', EPS, ' at T,H = ', T, H EPS = EPS*DABS(H03) C write (*,'(1x,E20.12)') Q, EPS IF (EPS.LE.DELTA) GO TO 111 DO 112 I=0,N1 DO 112 J=0,N2 112 Y(I,J) = Y0(I,J) T = T0 GO TO 113 C here we insert the time step measurement 111 if (SYS_ID.EQ.1) then N_TST = N_TST + 1 T_STEP(N_TST) = H endif IF (.NOT.LAST) GO TO 114 IF (FIRST) STEP = H RETURN 114 IF (EPS.EQ.0.0) GO TO 101 113 H = (DELTA/EPS)**0.2D00*0.8D00*H C write(*,*) 'H changed to ', H STEP = H FIRST = .FALSE. GO TO 102 END C-----------------------------------------------------------------------C C Prepiseme do podoby odpovidajici prave strane rovnic C hodnoty prave strany se ukladaji do pole GG SUBROUTINE RIGHT_HAND_SIDE (T, Phi, GG) DOUBLE PRECISION Phi(-1:401,-1:401), GG(0:400,0:400), T, / Kn(-1:401,-1:401), DT, F, H1, H2, ZETA, THETA_0, E, / R, EPS, ALPHA, vv1, vv2, vv3, Psi, T0_1, T0_2, T0_3, a, b, c DOUBLE PRECISION QQ(0:401,0:401), v1(0:401,0:401),v2(0:401,0:401), / NN1(0:401,0:401), NN2(0:401,0:401), NN3(0:401,0:401), / NU1(0:401,0:401), NU2(0:401,0:401), NU3(0:401,0:401), / TH1(0:401,0:401), TH2(0:401,0:401), TH3(0:401,0:401) INTEGER WRITE_ID COMMON /NUMER/ H1, H2, DT, N1, N2, EPS /ID/ WRITE_ID / /PHYS/ ALPHA, F, ZETA, THETA_0, M_FOLD, R / /ANIS/ a, b, c C functions of anisotropy C the anisotropy rotated 45 deg on Sep 4, 2006 c Psi (vv1,vv2,vv3) = DSQRT((1.D00+a)*vv1**2 + b*vv2**2 + c*vv3**2) c / +DSQRT(a*vv1**2 + (1.D00+b)*vv2**2 + c*vv3**2) c / +DSQRT(a*vv1**2 + b*vv2**2 + (1.D00+c)*vv3**2) c T0_1(vv1,vv2,vv3) = Psi (vv1,vv2,vv3)*vv1*( c / (1.D00+a)/DSQRT((1.D00+a)*vv1**2 + b*vv2**2 + c*vv3**2+E) c / +a/DSQRT(a*vv1**2 + (1.D00+b)*vv2**2 + c*vv3**2+E) c / +a/DSQRT(a*vv1**2 + b*vv2**2 + (1.D00+c)*vv3**2+E)) c T0_2(vv1,vv2,vv3) = Psi (vv1,vv2,vv3)*vv2*( c / b/DSQRT((1.D00+a)*vv1**2 + b*vv2**2 + c*vv3**2+E) c / +(1.D00+b)/DSQRT(a*vv1**2 + (1.D00+b)*vv2**2 + c*vv3**2+E) c / +b/DSQRT(a*vv1**2 + b*vv2**2 + (1.D00+c)*vv3**2+E)) c T0_3(vv1,vv2,vv3) = Psi (vv1,vv2,vv3)*vv3*( c / c/DSQRT((1.D00+a)*vv1**2 + b*vv2**2 + c*vv3**2+E) c / +c/DSQRT(a*vv1**2 + (1.D00+b)*vv2**2 + c*vv3**2+E) c / +(1.D00+c)/DSQRT(a*vv1**2 + b*vv2**2 + (1.D00+c)*vv3**2+E)) Psi (vv1,vv2,vv3) = / DSQRT((1.D00+a)*(vv1+vv2)**2 + b*(vv1-vv2)**2 + c*vv3**2) / +DSQRT(a*(vv1+vv2)**2 + (1.D00+b)*(vv1-vv2)**2 + c*vv3**2) / +DSQRT(a*(vv1+vv2)**2 + b*(vv1-vv2)**2 + (1.D00+c)*vv3**2) T0_1(vv1,vv2,vv3) = Psi (vv1,vv2,vv3)*( / ((1.D00+a)*(vv1+vv2)+b*(vv1-vv2)) / /DSQRT((1.D00+a)*(vv1+vv2)**2 + b*(vv1-vv2)**2 + c*vv3**2+E) / +(a*(vv1+vv2)+(1.D00+b)*(vv1-vv2)) / /DSQRT(a*(vv1+vv2)**2 + (1.D00+b)*(vv1-vv2)**2 + c*vv3**2+E) / +(a*(vv1+vv2)+b*(vv1-vv2)) / /DSQRT(a*(vv1+vv2)**2 + b*(vv1-vv2)**2 + (1.D00+c)*vv3**2+E)) T0_2(vv1,vv2,vv3) = Psi (vv1,vv2,vv3)*( / ((1.D00+a)*(vv1+vv2)-b*(vv1-vv2)) / /DSQRT((1.D00+a)*(vv1+vv2)**2 + b*(vv1-vv2)**2 + c*vv3**2+E) / +(a*(vv1+vv2)-(1.D00+b)*(vv1-vv2)) / /DSQRT(a*(vv1+vv2)**2 + (1.D00+b)*(vv1-vv2)**2 + c*vv3**2+E) / +(a*(vv1+vv2)-b*(vv1-vv2)) / /DSQRT(a*(vv1+vv2)**2 + b*(vv1-vv2)**2 + (1.D00+c)*vv3**2+E)) T0_3(vv1,vv2,vv3) = Psi (vv1,vv2,vv3)*vv3*( / c/DSQRT((1.D00+a)*(vv1+vv2)**2 + b*(vv1-vv2)**2 + c*vv3**2+E) / +c/DSQRT(a*(vv1+vv2)**2 + (1.D00+b)*(vv1-vv2)**2 + c*vv3**2+E) / +(1.D00+c)/ / DSQRT(a*(vv1+vv2)**2+b*(vv1-vv2)**2+(1.D00+c)*vv3**2+E)) E = 1.D-9 C reflexion of Phi C edges and corners DO I=0,N1 Phi(I,-1) = Phi(I,0) Phi(I,N2+1) = Phi(I,N2) ENDDO DO J=0,N2 Phi(-1,J) = Phi(0,J) Phi(N1+1,J) = Phi(N1,J) ENDDO Phi(-1,-1) = Phi(0,0) Phi(N1+1,-1) = Phi(N1,0) Phi(-1,N2+1) = Phi(0,N2) Phi(N1+1,N2+1)= Phi(N1,N2) C preparing gradients and normals DO i1=0,N1+1 DO i2=0,N2+1 v1(i1,i2) = (Phi(i1,i2)-Phi(i1-1,i2))/h1 v2(i1,i2) = (Phi(i1,i2)-Phi(i1,i2-1))/h2 C QQ(i1,i2) = DSQRT(eps**2+v1(i1,i2)**2+v2(i1,i2)**2) QQ(i1,i2) = Psi(-v1(i1,i2),-v2(i1,i2),1.D00) NU1(i1,i2) = -v1(i1,i2)/QQ(i1,i2) NU2(i1,i2) = -v2(i1,i2)/QQ(i1,i2) NU3(i1,i2) = 1.D00/QQ(i1,i2) NN1(i1,i2) = T0_1(NU1(i1,i2),NU2(i1,i2),NU3(i1,i2)) NN2(i1,i2) = T0_2(NU1(i1,i2),NU2(i1,i2),NU3(i1,i2)) NN3(i1,i2) = T0_3(NU1(i1,i2),NU2(i1,i2),NU3(i1,i2)) ENDDO ENDDO C computing curvature DO i1=0,N1 DO i2=0,N2 Kn(i1,i2) = / -(NN1(i1+1,i2)-NN1(i1,i2))/h1-(NN2(i1,i2+1)-NN2(i1,i2))/h2 ENDDO ENDDO C reflexion of curvature DO I=0,N1 Kn(I,-1) = Kn(I,0) Kn(I,N2+1) = Kn(I,N2) ENDDO DO J=0,N2 Kn(-1,J) = Kn(0,J) Kn(N1+1,J) = Kn(N1,J) ENDDO Kn(-1,-1) = Kn(0,0) Kn(N1+1,-1) = Kn(N1,0) Kn(-1,N2+1) = Kn(0,N2) Kn(N1+1,N2+1) = Kn(N1,N2) C preparing gradients C adding force according to Eberhard and previous codes in 1D DO i1=0,N1+1 DO i2=0,N2+1 v1(i1,i2) = (Kn(i1,i2)-Kn(i1-1,i2))/h1 / +F*(1.D00/Phi(i1,i2)-1.D00/Phi(i1-1,i2))/h1 v2(i1,i2) = (Kn(i1,i2)-Kn(i1,i2-1))/h2 / +F*(1.D00/Phi(i1,i2)-1.D00/Phi(i1,i2-1))/h2 TH1(i1,i2) = T0_1(v1(i1,i2),v2(i1,i2),0.0D00) TH2(i1,i2) = T0_2(v1(i1,i2),v2(i1,i2),0.0D00) TH3(i1,i2) = T0_3(v1(i1,i2),v2(i1,i2),0.0D00) ENDDO ENDDO C time evolution DO I=0,N1 DO J=0,N2 GG(I,J) = -(QQ(I+1,J)*(TH1(I+1,J)- / (TH1(I+1,J)*NU1(I+1,J)+TH2(I+1,J)*NU2(I+1,J) / +TH3(I+1,J)*NU3(I+1,J))*NN1(I+1,J))- / QQ(I,J)*(TH1(I,J)- / (TH1(I,J)*NU1(I,J)+TH2(I,J)*NU2(I,J) / +TH3(I,J)*NU3(I,J))*NN1(I,J)))/H1 / -(QQ(I,J+1)*(TH2(I,J+1)- / (TH1(I,J+1)*NU1(I,J+1)+TH2(I,J+1)*NU2(I,J+1) / +TH3(I,J+1)*NU3(I,J+1))*NN2(I,J+1))- / QQ(I,J)*(TH2(I,J)- / (TH1(I,J)*NU1(I,J)+TH2(I,J)*NU2(I,J) / +TH3(I,J)*NU3(I,J))*NN2(I,J)))/H2 ENDDO ENDDO RETURN END C-----------------------------------------------------------------------C SUBROUTINE SAVE_FIELD (Phi, IT) DOUBLE PRECISION Phi(-1:401,-1:401), H1, H2, DT, EPS, / Kn(-1:401,-1:401), GRADX(0:400,0:400), GRADY(0:400,0:400), / GR_NORM(0:400,0:400), QQ(0:400,0:400), QQ_m(0:400,0:400), / Dx_0, Dy_0, Dx_p, Dy_p, Dx_m, Dy_m, Dxx, Dyy, Dxy, NABLA_p, / NABLA_m DOUBLE PRECISION F, DF, DDF, R, K, CURV(0:400,0:400), S, CURV_AUX CHARACTER*12 FNAME_NUM, FNAME_EXA CHARACTER*4 TLEV_INDX COMMON /NUMER/ H1, H2, DT, N1, N2, EPS C----------------------------------------------------------------------C encode (4,'(i4.4)', TLEV_INDX) IT FNAME_NUM = 'f_'//TLEV_INDX//'.dat' WRITE(*,*) 'Writing computed temperature field at time level = ', / IT,' to the file ', FNAME_NUM OPEN (77,FILE=FNAME_NUM,STATUS = 'UNKNOWN') 10 FORMAT (17H# Time level T = , F8.5, 2I10) 12 FORMAT (2x,20F5.1) WRITE (77,10) IT*DT C prepiseme vystup do tvaru s (I*H1), Phi(1,I), Phi(2,I), treti souradnice DO I=0,N1 DO J=0,N2 WRITE (77,'(1X,2F15.5,2F25.5)') I*H1, J*H2, Phi(I,J), Kn(I,J) ENDDO WRITE (77,*) ENDDO CLOSE (77) RETURN END C-----------------------------------------------------------------------C C Hills shape C zadani pocatecni podminky - bude zadani Phi jako vektorove fce v uzlech SUBROUTINE INITIAL_LEVEL (Phi) DOUBLE PRECISION Phi(-1:401,-1:401), H1, H2, DT, R, X, Y, NORM, / F, ZETA, THETA_0, EPS, ALPHA, SCALE, Pi, K, R_ax, / THETA COMMON /NUMER/ H1, H2, DT, N1, N2, EPS / /PHYS/ ALPHA, F, ZETA, THETA_0, M_FOLD, R / /PATTERN/ SCALE /AX/ R_ax NORM(X,Y) = DSQRT(X*X+Y*Y) Pi = 3.1415926536D00 K = 2.D00 C Signed distance function smoothed by cut-off C-----------------------------------------------------------------------C C initial convex perturbation of the phase field C C Curve definition C C C C x_1 = R*cos(t), C C x_2 = R*sin(t). C C C C domain described by the condition [x_1-NUCL_1,x_2-NUCL_2] < curve C C H_crit_1, H_crit_2 have no meaning here C C sinusoidal perturbation I0 = N1/2 J0 = N2/2 DO I=-1,N1 DO J=-1,N2 THETA = DATAN2((J-J0)*H2,(I-I0)*H1+1.D-9) RR = NORM(DBLE(I-I0)*H1,DBLE(J-J0)*H2) C RR = NORM(DBLE(I-I0)*H1,0.D00) C RR = NORM(DBLE(I)*H1,0.D00)+SCALE*DCOS(2.*Pi*J/N2) C Phi(I,J) = 1.D00 + ZETA*DSIN(ALPHA*(RR-R))* C / DEXP(-ALPHA**2*(RR-R)**2)*(THETA_0+DSIGN(1.D00,RR-R)) C theta_0 is the mountain frequency, localized in a circle Phi(I,J) = 1.D00 + / K*DSIN(2.D00*PI*THETA_0*I/DBLE(N1))**M_FOLD* / DSIN(2.D00*PI*THETA_0*J/DBLE(N2))**M_FOLD* / (1.D00-DTANH(ALPHA*(RR-R)))/2.D00 ENDDO ENDDO C-----------------------------------------------------------------------C C edges and corners c DO I=0,N1 c Phi(I,-1) = Phi(I,1) c Phi(I,N2+1) = Phi(I,N2-1) c ENDDO c DO J=0,N2 c Phi(-1,J) = Phi(N1-1,J) c Phi(N1+1,J) = Phi(1,J) c ENDDO c Phi(-1,-1) = Phi(-1,1) c Phi(N1+1,-1) = Phi(N1+1,1) c Phi(-1,N2+1) = Phi(-1,N2-1) c Phi(N1+1,N2+1)= Phi(N1+1,N2-1) RETURN END C-----------------------------------------------------------------------C C loading and saving data SUBROUTINE STORAGE (Phi) DOUBLE PRECISION Phi(-1:401,-1:401), H1, H2, DT, R, X, Y, NORM, / F, ZETA, THETA_0, EPS, ALPHA, SCALE, Pi, K, / STEP INTEGER SAVE_ID, LOAD_ID COMMON /NUMER/ H1, H2, DT, N1, N2, EPS /MERSN/ STEP / /PHYS/ ALPHA, F, ZETA, THETA_0, M_FOLD, R / /PATTERN/ SCALE /STORE/ SAVE_ID, LOAD_ID IF (LOAD_ID.EQ.1) THEN LOAD_ID = 0 OPEN (77,FILE='elast22.ini',STATUS = 'UNKNOWN') READ (77,'(1X,F30.15)') STEP DO I=0,N1 READ (77,'(1X,20F30.15)') (Phi(I,J), J=0,N2) ENDDO CLOSE (77) WRITE (*,*) 'Data loaded from elast22.ini' ELSE OPEN (77,FILE='elast22.ini',STATUS = 'UNKNOWN') WRITE (77,'(1X,F30.15)') STEP DO I=0,N1 WRITE (77,'(1X,20F30.15)') (Phi(I,J), J=0,N2) ENDDO CLOSE (77) WRITE (*,*) 'Data saved into elast22.ini' ENDIF RETURN END C-----------------------------------------------------------------------C