program perturbsm c The intent of this program is to create perturbations c of soil moisture based on soil moisture data from the c Eta,GFS and RUC models. The RUC soil moisture was c interpolated from the 5 levels to 4 layers and then the c residual value was added to be consistent with the other c soil moisture sources. Qmin was then subtracted in the c ruc lsm in the /phys dir. parameter(im=256,jm=250,lm=4) real soilm005_ruc(im,jm),soilm020_ruc(im,jm) &, soilm040_ruc(im,jm),soilm160_ruc(im,jm) &, soilm300_ruc(im,jm),sm010200_gfs(im,jm) &, sm000010_eta(im,jm),sm000010_gfs(im,jm) &, sm010040_eta(im,jm),sm010040_gfs(im,jm) &, sm040100_eta(im,jm),sm040100_gfs(im,jm) &, sm100200_eta(im,jm),sm100200_gfs(im,jm) &, sm000010_ruc(im,jm),sm010040_ruc(im,jm) &, sm040100_ruc(im,jm),sm100200_ruc(im,jm) &, sm_ruc(im,jm,lm),sm_eta(im,jm,lm),sm_gfs(im,jm,lm) &, qmin(im,jm),soiltype(im,jm) real perturb1(im,jm,lm),perturb2(im,jm,lm) &, perturb3(im,jm,lm),perturb4(im,jm,lm) &, perturb5(im,jm,lm),perturb6(im,jm,lm) &, perturb7(im,jm,lm),perturb8(im,jm,lm) &, perturb9(im,jm,lm),perturb10(im,jm,lm) &, xbarminusx2_ruc(im,jm,lm),xbarminusx2_eta(im,jm,lm) &, xbarminusx2_gfs(im,jm,lm),mean(im,jm,lm) &, sum(im,jm,lm) real q1(im,jm,lm),q2(im,jm,lm),q3(im,jm,lm) &, soilm005_rucorig(im,jm),soilm020_rucorig(im,jm) &, soilm040_rucorig(im,jm),soilm160_rucorig(im,jm) &, soilm300_rucorig(im,jm) character*80 yearin integer year open(1,file='soilm005_ruc') open(2,file='soilm020_ruc') open(3,file='soilm040_ruc') open(4,file='soilm160_ruc') open(5,file='soilm300_ruc') do j=1,jm do i=1,im read(1,*)soilm005_ruc(i,j) enddo enddo do j=1,jm do i=1,im read(2,*)soilm020_ruc(i,j) enddo enddo do j=1,jm do i=1,im read(3,*)soilm040_ruc(i,j) enddo enddo do j=1,jm do i=1,im read(4,*)soilm160_ruc(i,j) enddo enddo do j=1,jm do i=1,im read(5,*)soilm300_ruc(i,j) enddo enddo close(1) close(2) close(3) close(4) close(5) c Read in RUC soil type open(1,file='rucsoiltype.txt') do j=1,jm do i=1,im read(1,*)soiltype(i,j) enddo enddo close(1) c Assign qmin to each grid point c Note that the #5,7,10 and 15 soil types don't c exist for my MS domain. I included it below, but c assumed it to be zero. do j=1,jm do i=1,im if(soiltype(i,j).eq.1)then qmin(i,j)=0.045 elseif(soiltype(i,j).eq.2)then qmin(i,j)=0.057 elseif(soiltype(i,j).eq.3)then qmin(i,j)=0.065 elseif(soiltype(i,j).eq.4)then qmin(i,j)=0.067 elseif(soiltype(i,j).eq.5)then qmin(i,j)=0. elseif(soiltype(i,j).eq.6)then qmin(i,j)=0.078 elseif(soiltype(i,j).eq.7)then qmin(i,j)=0. elseif(soiltype(i,j).eq.8)then qmin(i,j)=0.089 elseif(soiltype(i,j).eq.9)then qmin(i,j)=0.095 elseif(soiltype(i,j).eq.10)then qmin(i,j)=0. elseif(soiltype(i,j).eq.11)then qmin(i,j)=0.07 elseif(soiltype(i,j).eq.12)then qmin(i,j)=0.068 elseif(soiltype(i,j).eq.13)then qmin(i,j)=0.078 elseif(soiltype(i,j).eq.14)then qmin(i,j)=0. elseif(soiltype(i,j).eq.15)then qmin(i,j)=0. elseif(soiltype(i,j).eq.16)then qmin(i,j)=0.0649 endif enddo enddo print*,'qmin is ',qmin(128,125) c Add qmin(soiltype) to the RUC soil moisture values do j=1,jm do i=1,im soilm005_rucorig(i,j)=soilm005_ruc(i,j) soilm020_rucorig(i,j)=soilm020_ruc(i,j) soilm040_rucorig(i,j)=soilm040_ruc(i,j) soilm160_rucorig(i,j)=soilm160_ruc(i,j) soilm300_rucorig(i,j)=soilm300_ruc(i,j) soilm005_ruc(i,j)=soilm005_ruc(i,j)+qmin(i,j) soilm020_ruc(i,j)=soilm020_ruc(i,j)+qmin(i,j) soilm040_ruc(i,j)=soilm040_ruc(i,j)+qmin(i,j) soilm160_ruc(i,j)=soilm160_ruc(i,j)+qmin(i,j) soilm300_ruc(i,j)=soilm300_ruc(i,j)+qmin(i,j) enddo enddo c Interpolate RUC LSM levels to NOAH LSM levels c NOAH LSM RUC LSM c c ------------ 5cm 5cm ------------ c ------------ 25cm 20cm ------------ c 40cm ------------ c ------------ 70cm c c c ------------ 150cm c 160cm------------ c c c c c 300cm------------ c c 0-10cm use 5cm RUC c 10-40cm use y=y0 + alpha(y1-y0) c alpha=(x-x0/x1-x0) c x=soil level interested in c x0=soil level above x c x1=soil level below x c y=soilm020 + 0.25(soilm040 - soilm20) c 40-100cm c y=soilm040 + 0.25(soilm160-soilm40) c 100-200cm c y=soilm40 + (11./12.)(soilm160-soilm40) do j=1,jm do i=1,im sm000010_ruc(i,j)=soilm005_ruc(i,j) sm010040_ruc(i,j)=soilm020_ruc(i,j)+ 1 (0.25*(soilm040_ruc(i,j)-soilm020_ruc(i,j))) sm040100_ruc(i,j)=soilm040_ruc(i,j)+ 1 (0.25*(soilm160_ruc(i,j)-soilm040_ruc(i,j))) sm100200_ruc(i,j)=soilm040_ruc(i,j)+ 1 ((11./12.)*(soilm160_ruc(i,j)-soilm040_ruc(i,j))) enddo enddo print*,'ruc sm and ruc sm+qmin for 5 cm:' &, soilm005_rucorig(128,125) &, soilm005_ruc(128,125) print*,'ruc sm and ruc sm+qmin for 20 cm:' &, soilm020_rucorig(128,125) &, soilm020_ruc(128,125) print*,'ruc sm and ruc sm+qmin for 40 cm:' &, soilm040_rucorig(128,125) &, soilm040_ruc(128,125) print*,'ruc sm and ruc sm+qmin for 160 cm:' &, soilm160_rucorig(128,125) &, soilm160_ruc(128,125) print*,'ruc sm and ruc sm+qmin for 300 cm:' &, soilm300_rucorig(128,125) &, soilm300_ruc(128,125) c print*,'orig rucsm,rucsm+qmin,rucsminterp' c &, soilm005_rucorig(128,125) c &, soilm005_ruc(128,125),sm000010_ruc(128,125) print*,'Interpolated ruc sm: ' &, sm000010_ruc(128,125),sm010040_ruc(128,125) &, sm040100_ruc(128,125),sm100200_ruc(128,125) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Find year. Prior to June 2005 the GFS model had only two soil c layers (0-10cm and 10-200cm). After that time c the 4 layer soil fields were developed. If using data prior c to June 2005 this program will linearly interpolate c the GFS soil moisture to the NOAH LSM. The last layer soil c moisture will be assigned the previous layer's soil moisture cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc call getarg(1,yearin) read(yearin,*)year print*,'Year is: ',year c Read in GFS and Eta soil moisture data from wrf_real files if(year.ne.2005)then open(1,file='sm000010_eta') open(2,file='sm010040_eta') open(3,file='sm040100_eta') open(4,file='sm100200_eta') open(10,file='sm000010_gfs') open(40,file='sm010200_gfs') do j=1,jm do i=1,im read(10,*)sm000010_gfs(i,j) enddo enddo do j=1,jm do i=1,im read(40,*)sm010200_gfs(i,j) enddo enddo else open(1,file='sm000010_eta') open(2,file='sm010040_eta') open(3,file='sm040100_eta') open(4,file='sm100200_eta') open(10,file='sm000010_gfs') open(20,file='sm010040_gfs') open(30,file='sm040100_gfs') open(40,file='sm100200_gfs') do j=1,jm do i=1,im read(10,*)sm000010_gfs(i,j) enddo enddo do j=1,jm do i=1,im read(20,*)sm010040_gfs(i,j) enddo enddo do j=1,jm do i=1,im read(30,*)sm040100_gfs(i,j) enddo enddo do j=1,jm do i=1,im read(40,*)sm100200_gfs(i,j) enddo enddo endif do j=1,jm do i=1,im read(1,*)sm000010_eta(i,j) enddo enddo do j=1,jm do i=1,im read(2,*)sm010040_eta(i,j) enddo enddo do j=1,jm do i=1,im read(3,*)sm040100_eta(i,j) enddo enddo do j=1,jm do i=1,im read(4,*)sm100200_eta(i,j) enddo enddo if(year.ne.2005)then close(1) close(2) close(3) close(4) close(10) close(40) else close(1) close(2) close(3) close(4) close(10) close(20) close(30) close(40) endif print*,'Now the Eta: ' print*,'soil moisture for 0-10cm: ' &, sm000010_eta(128,125) print*,'soil moisture for 10-40cm: ' &, sm010040_eta(128,125) print*,'soil moisture for 40-100cm: ' &, sm040100_eta(128,125) print*,'soil moisture for 100-200cm: ' &, sm100200_eta(128,125) print*,'Now the GFS: ' if(year.ne.2005)then print*,'ONLY 2 SOIL LAYERS' print*,'soil moisture for 0-10cm: ' &, sm000010_gfs(128,125) print*,'soil moisture for 10-200cm: ' &, sm010200_gfs(128,125) else print*,'soil moisture for 0-10cm: ' &, sm000010_gfs(128,125) print*,'soil moisture for 10-40cm: ' &, sm010040_gfs(128,125) print*,'soil moisture for 40-100cm: ' &, sm040100_gfs(128,125) print*,'soil moisture for 100-200cm: ' &, sm100200_gfs(128,125) endif c NOAH LSM GFS c c --------------5cm 5cm------------- c c --------------25cm x c c c---------------70cm x c c 105cm---------- c c---------------150cm c y=y0+alpha(y1-y0) c alpha=(x-x0/x1-x0) c For 10-40cm,alpha=(25-5)/(105-5)=0.2 c For 40-100cm,alpha=(70-5)/(105-5)=0.65 if(year.ne.2005)then do j=1,jm do i=1,im sm010040_gfs(i,j)=sm000010(i,j)+ 1 0.2*(sm010200(i,j)-sm000010(i,j)) sm040100_gfs(i,j)=sm000010(i,j)+ 1 0.65*(sm010200(i,j)-sm000010(i,j)) sm100200_gfs(i,j)=sm040100_gfs(i,j) enddo enddo endif c combine one level data into multiple levels do j=1,jm do i=1,im sm_ruc(i,j,1)=sm000010_ruc(i,j) sm_ruc(i,j,2)=sm010040_ruc(i,j) sm_ruc(i,j,3)=sm040100_ruc(i,j) sm_ruc(i,j,4)=sm100200_ruc(i,j) sm_eta(i,j,1)=sm000010_eta(i,j) sm_eta(i,j,2)=sm010040_eta(i,j) sm_eta(i,j,3)=sm040100_eta(i,j) sm_eta(i,j,4)=sm100200_eta(i,j) sm_gfs(i,j,1)=sm000010_gfs(i,j) sm_gfs(i,j,2)=sm010040_gfs(i,j) sm_gfs(i,j,3)=sm040100_gfs(i,j) sm_gfs(i,j,4)=sm100200_gfs(i,j) enddo enddo c Calculate the soil moisture perturbations for each grid point c and one level at a time do l=1,lm do j=1,jm do i=1,im ccc First calculate mean mean(i,j,l)=(sm_ruc(i,j,l)+sm_eta(i,j,l) 1 +sm_gfs(i,j,l))/3. ccc RUC c perturb1(i,j,l)=0.0 perturb1(i,j,l)=sm_ruc(i,j,l) ccc ETA perturb2(i,j,l)=sm_eta(i,j,l) ccc GFS perturb3(i,j,l)=sm_gfs(i,j,l) ccc standard deviations xbarminusx2_ruc(i,j,l)=(mean(i,j,l)- 1 sm_ruc(i,j,l))**2. xbarminusx2_eta(i,j,l)=(mean(i,j,l)- 1 sm_eta(i,j,l))**2. xbarminusx2_gfs(i,j,l)=(mean(i,j,l)- 1 sm_gfs(i,j,l))**2. sum(i,j,l)=xbarminusx2_ruc(i,j,l)+xbarminusx2_eta(i,j,l) 1 + xbarminusx2_gfs(i,j,l) perturb4(i,j,l)= 1.*sqrt(sum(i,j,l)/2.)+mean(i,j,l) perturb5(i,j,l)= -1.*sqrt(sum(i,j,l)/2.)+mean(i,j,l) perturb6(i,j,l)= 2.*sqrt(sum(i,j,l)/2.)+mean(i,j,l) perturb7(i,j,l)= -2.*sqrt(sum(i,j,l)/2.)+mean(i,j,l) caligo c if(perturb7(i,j,l).lt.0)then c perturb7(i,j,l)=0. c endif ccc mean perturb8(i,j,l)=mean(i,j,l) ccc percentiles ccccc arrange values from least to greatest q1(i,j,l)=min(sm_ruc(i,j,l),sm_eta(i,j,l),sm_gfs(i,j,l)) if(q1(i,j,l).eq.sm_ruc(i,j,l))then q2(i,j,l)=min(sm_eta(i,j,l),sm_gfs(i,j,l)) elseif(q1(i,j,l).eq.sm_eta(i,j,l))then q2(i,j,l)=min(sm_ruc(i,j,l),sm_gfs(i,j,l)) elseif(q1(i,j,l).eq.sm_gfs(i,j,l))then q2(i,j,l)=min(sm_ruc(i,j,l),sm_eta(i,j,l)) endif q3(i,j,l)=max(sm_ruc(i,j,l),sm_eta(i,j,l),sm_gfs(i,j,l)) ccccc calculate p values (needed only once) pvalue1=(1-.5)/3. ! 0.1666 pvalue2=(2-.5)/3. ! 0.50 pvalue3=(3-.5)/3. ! 0.8333 ccccc Q(0.25) c f=n*(p - pvalue1),where n is # of obs,p=p-value c interested in and pvalue1 is one p-value lower f=3.*(0.25-pvalue1) perturb9(i,j,l)= (1.-f)*(q1(i,j,l)) + f*(q2(i,j,l)) ccccc Q(0.75) f=3.*(0.75-pvalue2) perturb10(i,j,l)= (1.-f)*(q2(i,j,l)) + f*(q3(i,j,l)) enddo enddo enddo do l=1,lm print*,'Mean,l : ',mean(128,125,l) print*,'25th percentile,l: ',perturb9(128,125,l) print*,'75th percentile,l: ',perturb10(128,125,l) print*,' ' print*,'q1,q2,q3,l: ',q1(128,125,l),q2(128,125,l) &,q3(128,125,l) print*,' ' print*,'1sigma,-1sigma,2sigma,-2sigma:' &,perturb4(128,125,l),perturb5(128,125,l),perturb6(128,125,l) &,perturb7(128,125,l) print*,' ' print*,'sbarminusx2(ruc,eta,gfs),sum: ' &,xbarminusx2_ruc(128,125,l),xbarminusx2_eta(128,125,l) &,xbarminusx2_gfs(128,125,l) enddo c Write out perturbations in special format cc Peturbation 1 do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(16,*)perturb1(i,j,1),';' else write(16,*)perturb1(i,j,1),',' endif enddo enddo c here you would put more variables.... c c c do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(16,*)perturb1(i,j,2),';' else if(i.eq.1.and.j.eq.1)then write(16,*)'SM010040=' endif write(16,*)perturb1(i,j,2),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(16,*)perturb1(i,j,3),';' else if(i.eq.1.and.j.eq.1)then write(16,*)'SM040100=' endif write(16,*)perturb1(i,j,3),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(16,*)perturb1(i,j,4),';' else if(i.eq.1.and.j.eq.1)then write(16,*)'SM100200=' endif write(16,*)perturb1(i,j,4),',' endif enddo enddo c go to 100 c go to 111 C Write out perturbation 2 to a file do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(17,*)perturb2(i,j,1),';' else write(17,*)perturb2(i,j,1),',' endif enddo enddo c here you would put more variables.... c c c do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(17,*)perturb2(i,j,2),';' else if(i.eq.1.and.j.eq.1)then write(17,*)'SM010040=' endif write(17,*)perturb2(i,j,2),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(17,*)perturb2(i,j,3),';' else if(i.eq.1.and.j.eq.1)then write(17,*)'SM040100=' endif write(17,*)perturb2(i,j,3),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(17,*)perturb2(i,j,4),';' else if(i.eq.1.and.j.eq.1)then write(17,*)'SM100200=' endif write(17,*)perturb2(i,j,4),',' endif enddo enddo C Peturbation 3 do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(18,*)perturb3(i,j,1),';' else write(18,*)perturb3(i,j,1),',' endif enddo enddo c here you would put more variables.... c c c do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(18,*)perturb3(i,j,2),';' else if(i.eq.1.and.j.eq.1)then write(18,*)'SM010040=' endif write(18,*)perturb3(i,j,2),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(18,*)perturb3(i,j,3),';' else if(i.eq.1.and.j.eq.1)then write(18,*)'SM040100=' endif write(18,*)perturb3(i,j,3),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(18,*)perturb3(i,j,4),';' else if(i.eq.1.and.j.eq.1)then write(18,*)'SM100200=' endif write(18,*)perturb3(i,j,4),',' endif enddo enddo 111 continue C Peturbation 4 do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(19,*)perturb4(i,j,1),';' else write(19,*)perturb4(i,j,1),',' endif enddo enddo c here you would put more variables.... c c c do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(19,*)perturb4(i,j,2),';' else if(i.eq.1.and.j.eq.1)then write(19,*)'SM010040=' endif write(19,*)perturb4(i,j,2),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(19,*)perturb4(i,j,3),';' else if(i.eq.1.and.j.eq.1)then write(19,*)'SM040100=' endif write(19,*)perturb4(i,j,3),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(19,*)perturb4(i,j,4),';' else if(i.eq.1.and.j.eq.1)then write(19,*)'SM100200=' endif write(19,*)perturb4(i,j,4),',' endif enddo enddo C Peturbation 5 c Write out RUC,ETA and GFS values in special format do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(20,*)perturb5(i,j,1),';' else write(20,*)perturb5(i,j,1),',' endif enddo enddo c here you would put more variables.... c c c do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(20,*)perturb5(i,j,2),';' else if(i.eq.1.and.j.eq.1)then write(20,*)'SM010040=' endif write(20,*)perturb5(i,j,2),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(20,*)perturb5(i,j,3),';' else if(i.eq.1.and.j.eq.1)then write(20,*)'SM040100=' endif write(20,*)perturb5(i,j,3),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(20,*)perturb5(i,j,4),';' else if(i.eq.1.and.j.eq.1)then write(20,*)'SM100200=' endif write(20,*)perturb5(i,j,4),',' endif enddo enddo C Peturbation 6 do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(21,*)perturb6(i,j,1),';' else write(21,*)perturb6(i,j,1),',' endif enddo enddo c here you would put more variables.... c c c do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(21,*)perturb6(i,j,2),';' else if(i.eq.1.and.j.eq.1)then write(21,*)'SM010040=' endif write(21,*)perturb6(i,j,2),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(21,*)perturb6(i,j,3),';' else if(i.eq.1.and.j.eq.1)then write(21,*)'SM040100=' endif write(21,*)perturb6(i,j,3),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(21,*)perturb6(i,j,4),';' else if(i.eq.1.and.j.eq.1)then write(21,*)'SM100200=' endif write(21,*)perturb6(i,j,4),',' endif enddo enddo C Peturbation 7 do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(22,*)perturb7(i,j,1),';' else write(22,*)perturb7(i,j,1),',' endif enddo enddo c here you would put more variables.... c c c do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(22,*)perturb7(i,j,2),';' else if(i.eq.1.and.j.eq.1)then write(22,*)'SM010040=' endif write(22,*)perturb7(i,j,2),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(22,*)perturb7(i,j,3),';' else if(i.eq.1.and.j.eq.1)then write(22,*)'SM040100=' endif write(22,*)perturb7(i,j,3),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(22,*)perturb7(i,j,4),';' else if(i.eq.1.and.j.eq.1)then write(22,*)'SM100200=' endif write(22,*)perturb7(i,j,4),',' endif enddo enddo C Peturbation 8 do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(23,*)perturb8(i,j,1),';' else write(23,*)perturb8(i,j,1),',' endif enddo enddo c here you would put more variables.... c c c do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(23,*)perturb8(i,j,2),';' else if(i.eq.1.and.j.eq.1)then write(23,*)'SM010040=' endif write(23,*)perturb8(i,j,2),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(23,*)perturb8(i,j,3),';' else if(i.eq.1.and.j.eq.1)then write(23,*)'SM040100=' endif write(23,*)perturb8(i,j,3),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(23,*)perturb8(i,j,4),';' else if(i.eq.1.and.j.eq.1)then write(23,*)'SM100200=' endif write(23,*)perturb8(i,j,4),',' endif enddo enddo C Peturbation 9 do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(24,*)perturb9(i,j,1),';' else write(24,*)perturb9(i,j,1),',' endif enddo enddo c here you would put more variables.... c c c do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(24,*)perturb9(i,j,2),';' else if(i.eq.1.and.j.eq.1)then write(24,*)'SM010040=' endif write(24,*)perturb9(i,j,2),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(24,*)perturb9(i,j,3),';' else if(i.eq.1.and.j.eq.1)then write(24,*)'SM040100=' endif write(24,*)perturb9(i,j,3),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(24,*)perturb9(i,j,4),';' else if(i.eq.1.and.j.eq.1)then write(24,*)'SM100200=' endif write(24,*)perturb9(i,j,4),',' endif enddo enddo C Peturbation 10 do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(25,*)perturb10(i,j,1),';' else write(25,*)perturb10(i,j,1),',' endif enddo enddo c here you would put more variables.... c c c do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(25,*)perturb10(i,j,2),';' else if(i.eq.1.and.j.eq.1)then write(25,*)'SM010040=' endif write(25,*)perturb10(i,j,2),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(25,*)perturb10(i,j,3),';' else if(i.eq.1.and.j.eq.1)then write(25,*)'SM040100=' endif write(25,*)perturb10(i,j,3),',' endif enddo enddo do j=1,jm do i=1,im if ( i.eq.im .and. j.eq.jm) then write(25,*)perturb10(i,j,4),';' else if(i.eq.1.and.j.eq.1)then write(25,*)'SM100200=' endif write(25,*)perturb10(i,j,4),',' endif enddo enddo c write(16,*)'}' 100 continue end