SLiMv2.1: SLiMv2.1/derivative.f90 Source File

00001       subroutine initialize(ider,diff,dt) 
00002 ! this sets up the derivative operator (in F-space) ider.
00003 ! 2 pi -> dx * nx  
00004       implicit none
00005       include 'param.f90'   
00006       integer*4 k1,j1,i,j,k,offk
00007       complex ider(nx/2+1,ny,nz_max,3)
00008       complex diff(nx/2+1,ny,nz_max)
00009       real dt,fac,z,zcorr
00010       do i=1,nx/2+1
00011          do k=0,nz
00012            z=6.959e10-(1000.-k)*dz
00013            zcorr=z/6.959e10 
00014            ider(i,:,k,1)= (i-1.)*(0.,1.)/(dx*nx*zcorr/(2.*pi))
00015          enddo  
00016          !stop
00017          fac=1
00018          if ((i .gt. 15) .and. (i .lt. 40)) fac=1-(i-15)/35.  
00019 !next is for f-mode plage 27.may.2008
00020          diff(i,:,:)=diff(i,:,:)-abs(ider(i,:,:,1))**2.19*6.17314e13/nu/2.*1.3
00021 !         diff(i,:,:)=diff(i,:,:)-abs(ider(i,:,:,1))**2.19*6.17314e13/nu/2.*fac/1.5 !for p1
00022       enddo
00023       do j=1,ny
00024          j1=j
00025          if (j .gt. (ny+1)/2) j1=j-ny 
00026          ider(:,j,:,2)= (j1-1)*(0.,1.)/(dy*ny/(2.*pi))
00027          fac=1
00028          if ((j .gt. 15) .and. (j .lt. 40)) fac=1-(j-15)/35.
00029          diff(:,j,:)=diff(:,j,:)-abs(ider(:,j,:,2))**2.19*6.17314e13/nu/2.*1.3
00030       enddo   
00031 
00032         do k=1,nz          
00033            if(kk(k) .gt. nz_ghostglobal-80) then  
00034               offk=float(kk(k)-nz_ghostglobal+80)
00035          !     diff(:,:,k)=diff(:,:,k)+diff(25,1,k)*exp((offk-16)/3.)  !/4. !54
00036               diff(:,:,k)=diff(:,:,k)-2.5e-11*(exp((offk-30)/10.))   !(offk-27)/20.)
00037          endif
00038         enddo   
00039 
00040         do k=1,nz
00041            if(kk(k) .lt. 200) then        
00042                          offk=200-kk(k)
00043 !                         diff(:,:,k)=diff(:,:,k)+diff(25,1,200)*exp(float(offk-64)/50.)   !16
00044                          diff(:,:,k)=diff(:,:,k)-1./240./nu*exp((offk-80)/25.)    !-50   /25.
00045            endif              
00046         enddo   
00047 
00048 ! semi explicit scheme
00049       diff=(exp(diff(:,:,:)*nu*2*dt)-1.)/(2.*dt)
00050       end subroutine initialize
00051 !=====================================================================
00052       function centre(datain)
00053       implicit none
00054       include 'param.f90' 
00055       real datain(nx,ny,nz_max),centre(nx,ny,nz_max)
00056       include 'terms3d.f90' 
00057       if(step .eq. 1) then
00058          centre(:,:,1:nz-1)=(datain(:,:,1:nz-1)+datain(:,:,2:nz))/2.
00059          centre(:,:,nz)=(3*datain(:,:,nz-1)-datain(:,:,nz-2))/2.
00060       endif
00061       if(step .eq. 2) then
00062          centre(:,:,2:nz)=(datain(:,:,1:nz-1)+datain(:,:,2:nz))/2.
00063          centre(:,:,1)=(2*datain(:,:,2)-datain(:,:,2))/2.
00064       endif
00065       return
00066       end function centre
00067 !=====================================================================
00068       function centrebc(datain)
00069       implicit none
00070       include 'param.f90' 
00071       real datain(nx,ny,nz_max),centrebc(nx,ny)
00072       include 'terms3d.f90' 
00073       if(step .eq. 1) then
00074          centrebc(:,:)=(datain(:,:,nz-1)+datain(:,:,nz))/2.
00075       endif
00076       if(step .eq. 2) then
00077          centrebc(:,:)=(datain(:,:,nz-1)+datain(:,:,nz))/2.
00078       endif
00079       return
00080       end function centrebc
00081 !====================================================================
00082       function centref(datain)
00083       implicit none
00084       include 'param.f90' 
00085       complex datain(nx/2+1,ny,nz_max),centref(nx/2+1,ny,nz_max)
00086       include 'terms3d.f90' 
00087       if(step .eq. 1) then
00088          centref(:,:,1:nz-2)=(datain(:,:,1:nz-2)+datain(:,:,2:nz-1))/2.
00089          centref(:,:,nz-1)=(3.*datain(:,:,nz-2)-datain(:,:,nz-1))/2.
00090       endif
00091       if(step .eq. 2) then
00092          centref(:,:,2:nz)=(datain(:,:,1:nz-1)+datain(:,:,2:nz))/2.
00093          centref(:,:,1)=(3.*datain(:,:,2)-datain(:,:,3))/2.
00094          centref(:,:,1)=0
00095       endif
00096       return
00097       end function centref
00098 !====================================================================
00099       function ddxc(datain)
00100       implicit none
00101       include 'param.f90' 
00102        complex datain(nx/2+1,ny,nz_max),ddxc(nx/2+1,ny,nz_max),temp(nx/2+1,ny,nz_max)
00103       include 'terms3d.f90' 
00104       temp=datain*ider(:,:,:,1)
00105       if(step .eq. 1) then
00106          ddxc(:,:,1:nz-1)=(temp(:,:,1:nz-1)+temp(:,:,2:nz))/2.
00107          ddxc(:,:,nz)=(temp(:,:,nz-1)+temp(:,:,nz))/2.
00108       endif
00109       if(step .eq. 2) then
00110          ddxc(:,:,2:nz)=(temp(:,:,1:nz-1)+temp(:,:,2:nz))/2.
00111          ddxc(:,:,1)=(temp(:,:,1))*0
00112       endif
00113       return
00114       end function ddxc
00115 !=====================================================================
00116       function ddxe(datain)
00117       implicit none
00118       include 'param.f90' 
00119        complex datain(nx/2+1,ny,nz_max),ddxe(nx/2+1,ny,nz_max)
00120       include 'terms3d.f90' 
00121       ddxe=datain*ider(:,:,:,1)
00122       return
00123       end function ddxe
00124 !=====================================================================
00125       function ddyc(datain)
00126       implicit none
00127       include 'param.f90'
00128       complex datain(nx/2+1,ny,nz_max),ddyc(nx/2+1,ny,nz_max),temp(nx/2+1,ny,nz_max)
00129       include 'terms3d.f90' 
00130       temp=datain*ider(:,:,:,2)
00131       if(step .eq. 1) then
00132          ddyc(:,:,1:nz-2)=(temp(:,:,1:nz-2)+temp(:,:,2:nz-1))/2.
00133          ddyc(:,:,nz-1)=(temp(:,:,nz-2)+temp(:,:,nz-1))/2
00134       endif
00135       if(step .eq. 2) then
00136          ddyc(:,:,2:nz)=(temp(:,:,1:nz-1)+temp(:,:,2:nz))/2.
00137          ddyc(:,:,1)=(temp(:,:,1))*0
00138       endif
00139       return
00140       end function ddyc
00141 !=====================================================================
00142       function ddye(datain)
00143       implicit none
00144       include 'param.f90'
00145       complex datain(nx/2+1,ny,nz_max),ddye(nx/2+1,ny,nz_max)
00146       include 'terms3d.f90' 
00147       ddye=datain*ider(:,:,:,2)
00148       return
00149       end function ddye
00150 
00151 !=====================================================================
00152       function ddzc(datain)
00153 !differentiation for result at cell centres
00154         implicit none
00155         include 'param.f90'
00156         include 'terms3d.f90' 
00157         real datain(nx,ny,nz_max),ddzc(nx,ny,nz_max),ddztemp(nx,ny,nz_max),true
00158         !LAX-Wendroff (2-step)
00159         if(step .eq.1) then
00160            ddzc(:,:,1:nz-1)=(datain(:,:,2:nz  )-datain(:,:,1:nz-1))/(dz) 
00161            ddzc(:,:,nz)=(datain(:,:,nz  )-datain(:,:,nz-1))/(dz)          
00162         else
00163            ddzc(:,:,2:nz)=(datain(:,:,2:nz  )-datain(:,:,1:nz-1))/(dz)
00164            ddzc(:,:,1)=(datain(:,:,2  )-datain(:,:,1))/(dz)*0
00165         endif
00166         return
00167       end function ddzc
00168 !=====================================================================
00169       function ddze(datain)
00170 !differeniation for result at cell edges
00171         implicit none
00172         include 'param.f90'
00173         include 'terms3d.f90' 
00174         real datain(nx,ny,nz_max),ddze(nx,ny,nz_max),ddztemp(nx,ny,nz_max),true
00175            ddze(:,:,2:nz-1)=(datain(:,:,3:nz  )-datain(:,:,1:nz-2))/(2.*dz) 
00176            ddze(:,:,nz)=(datain(:,:,nz)-datain(:,:,nz-1))/(dz)        
00177            ddze(:,:,1)=(datain(:,:,2  )-datain(:,:,1))/(dz)*0          
00178         return
00179       end function ddze
00180 !=====================================================================
00181       function ddzc_curl(datain)
00182 !differentiation for result at cell centres
00183         implicit none
00184         include 'param.f90'
00185         include 'terms3d.f90' 
00186         real datain(nx,ny,nz_max),ddzc_curl(nx,ny,nz_max),ddztemp(nx,ny,nz_max),true
00187         !LAX-Wendroff (2-step)
00188         if(step .eq.1) then
00189            ddzc_curl(:,:,1:nz-1)=(datain(:,:,2:nz  )-datain(:,:,1:nz-1))/(dz) 
00190            ddzc_curl(:,:,nz)=(datain(:,:,nz  )-datain(:,:,nz-1))/(dz)          
00191         else
00192            ddzc_curl(:,:,2:nz)=(datain(:,:,2:nz  )-datain(:,:,1:nz-1))/(dz)
00193            ddzc_curl(:,:,1)=(datain(:,:,2  )-datain(:,:,1))/(dz)*0
00194         endif
00195         ddzc_curl= ddzc_curl+datain/zz
00196         return
00197       end function ddzc_curl
00198 !=====================================================================
00199       function ddze_curl(datain)
00200 !differeniation for result at cell edges
00201         implicit none
00202         include 'param.f90'
00203         include 'terms3d.f90' 
00204         real datain(nx,ny,nz_max),ddze_curl(nx,ny,nz_max),ddztemp(nx,ny,nz_max),true
00205            ddze_curl(:,:,2:nz-1)=(datain(:,:,3:nz  )-datain(:,:,1:nz-2))/(2.*dz) 
00206            ddze_curl(:,:,nz)=(datain(:,:,nz)-datain(:,:,nz-1))/(dz)        
00207 !spherical
00208            ddze_curl= ddze_curl+datain/zz
00209            ddze_curl(:,:,1)=(datain(:,:,2  )-datain(:,:,1))/(dz)*0          
00210         return
00211       end function ddze_curl
00212 !=====================================================================
00213       function ddzc_div(datain)
00214 !differentiation for result at cell centres
00215         implicit none
00216         include 'param.f90'
00217         include 'terms3d.f90' 
00218         real datain(nx,ny,nz_max),ddzc_div(nx,ny,nz_max),ddztemp(nx,ny,nz_max),true
00219         !LAX-Wendroff (2-step)
00220         if(step .eq.1) then
00221            ddzc_div(:,:,1:nz-1)=(datain(:,:,2:nz  )-datain(:,:,1:nz-1))/(dz) 
00222            ddzc_div(:,:,nz)=(datain(:,:,nz  )-datain(:,:,nz-1))/(dz)          
00223         else
00224            ddzc_div(:,:,2:nz)=(datain(:,:,2:nz  )-datain(:,:,1:nz-1))/(dz)
00225            ddzc_div(:,:,1)=(datain(:,:,2  )-datain(:,:,1))/(dz)*0
00226         endif
00227 !spherical
00228            ddzc_div= ddzc_div+2* datain/zz
00229         return
00230       end function ddzc_div
00231 !=====================================================================
00232       function ddze_div(datain)
00233 !differeniation for result at cell edges
00234         implicit none
00235         include 'param.f90'
00236         include 'terms3d.f90' 
00237         real datain(nx,ny,nz_max),ddze_div(nx,ny,nz_max),ddztemp(nx,ny,nz_max),true
00238            ddze_div(:,:,2:nz-1)=(datain(:,:,3:nz  )-datain(:,:,1:nz-2))/(2.*dz) 
00239            ddze_div(:,:,nz)=(datain(:,:,nz)-datain(:,:,nz-1))/(dz)       
00240 !Spherical 
00241            ddze_div= ddze_div+2.*datain/zz
00242            ddze_div(:,:,1)=(datain(:,:,2  )-datain(:,:,1))/(dz)*0 
00243         return
00244       end function ddze_div
00245 !=====================================================================
00246       function d2fdz2(datain)
00247       implicit none
00248       include 'param.f90'
00249 !      include 'terms3d.f90'
00250       integer k,k1
00251       complex datain(nx/2+1,ny,nz_max),d2fdz2(nx/2+1,ny,nz_max)
00252       d2fdz2(:,:,2:nz-2)=(datain(:,:,3:nz-1)-2.*datain(:,:,2:nz-2)  &
00253                 +datain(:,:,1:nz-3))/(dz**2) 
00254       d2fdz2(:,:,1)=0
00255       d2fdz2(:,:,nz-1)=(datain(:,:,nz-2)-2.*datain(:,:,nz-1)  &
00256                 +datain(:,:,nz-2))/(dz**2) 
00257 !More diffuion at Boundaries
00258 !     do k=1,nz
00259 !       if(kk(k) .le. 400) d2fdz2(:,:,k)=d2fdz2(:,:,k)*exp((201.-kk(k))/200.)   !200 201 100
00260 !
00261 !       if(kk(k) .gt. nz_ghostglobal-40)   then
00262 !                 k1=kk(k)-(nz_ghostglobal-20)
00263 !                 d2fdz2(:,:,k)=d2fdz2(:,:,k)*(1+exp(-k1/10.)) 
00264 !         endif
00265 !     enddo     
00266      return
00267      end function d2fdz2