C FORZQ.f
C
C VIS-5D analysis function to compute Q vector Forcing 
C (in 1.E-18 m kg-1 s-1 )

      INTEGER FUNCTION USERFUNC( OUTGRID, OUTNL, OUTLOWLEV,
     *                           INGRID, NR, NC, NL, LOWLEV, MAXNL,
     *                           NVARS, NAMES,
     *                           DATE, TIME,
     *                           PROJECTION, PROJ_ARGS,
     *                           VERTICAL, VERT_ARGS )
C     ARGUMENTS:
	implicit none
      INTEGER NVARS
      INTEGER NR, NC, NL(NVARS), LOWLEV(NVARS), MAXNL
      REAL OUTGRID(NR,NC,MAXNL)
      INTEGER OUTNL, OUTLOWLEV
      REAL INGRID(NR,NC,MAXNL,NVARS)
      CHARACTER*8 NAMES(NVARS)
      INTEGER DATE, TIME
      INTEGER PROJECTION
      REAL PROJ_ARGS(*)
      INTEGER VERTICAL
      REAL VERT_ARGS(*)

C     LOCAL VARS:
      integer var, ir, ic, il,it,ig,iqv,iqu
      integer ilat,imap
      real lat,Pi,Rt,Per
      real q1(200,200,22),Q2(200,200,22)
      real forzq,p(22),r
      real g,cor,xincx,xincy,b20,b21,b22,b23,div,ds
      parameter (Pi=3.141596,Rt=6371227.,Per=86400.,g=9.807,r=287.04)

C     Specify number of levels in OUTGRID
      OUTNL = MAXNL
      OUTLOWLEV = 0

C     Find the variables
      iqu = -1
      iqv = -1
      it = -1
      ig = -1
      ilat = -1
      imap = -1      
      do var=1,nvars
         if (names(var) .eq. 'UQ') iqu = var
         if (names(var) .eq. 'VQ') iqv = var
         if (names(var) .eq. 'LATITCRS') ilat = var
         if (names(var) .eq. 'MAPFACCR') imap = var                
      end do
      
       if (ilat .eq. -1 .or. imap .eq. -1) then
         print *, "Couldn't find LATITCRS and MAPFACCR vars!"
         userfunc = 1
         return
      endif   
      
      if (projection.Eq.2) then
        ds=proj_args(6)*1000.        
      else if (projection.Eq.3) then  
        ds=proj_args(5)*1000. 
      else
        print *, "Uncomprensible projection data!"
        userfunc = 1
        return
      endif
              
c BUSCANDO SI YA TENEMOS UQ I VQ CALCULADOS
C     If uq, or vq not found, return error 1
      if (iqu .eq. -1 .or. iqv .eq. -1) then
	 print*,'Calculating UQ and VQ internally'
	else
	go to 143
      endif

      do var=1,nvars
         if (names(var) .eq. 'T') it = var
         if (names(var) .eq. 'H') ig = var
      end do
C   If T, or H not found, return error 1
      if (it .eq. -1 .or. ig .eq. -1) then
         print *, "Couldn't find T, and/or H variables!"
         userfunc = 1
         return
      endif


c Calculo de los niveles de presion

c INCREMENTO EN Y (= 2. * Brazo de malla)

        xincy=2.*ds

C  Bucle sobre todos los niveles de presion
      do 100 il=1,maxnl
	p(il)=1012.5*exp(vert_args(il)/-7.2)*100.
C CALCULO DEL VECTOR Q (Q1,Q2) (en unidades de 1E-10 SI):

         do ic=3,nc-2
           do ir=3,nr-2
C CALCULO DE LATITUD, CORIOLIS Y INCREMENTO EN X (2.* Brazo de malla)
		lat=ingrid(ir,ic,1,ilat)*Pi/180.
		cor=4.*Pi/Per*Sin(lat)
		xincx=2.*ds

c CALCULO DE Q1
	
          B20=(-ingrid(ir+1,ic+1,il,ig)+ingrid(ir-1,ic+1,il,ig)+
     & ingrid(ir+1,ic-1,il,ig)-ingrid(ir-1,ic-1,il,ig))/(xincx*xincy)
        B20=B20*(ingrid(ir,ic+1,il,it)-ingrid(ir,ic-1,il,it))/xincx
        B20=B20*(ingrid(ir,ic,1,imap)**3.)
          B21=(ingrid(ir,ic+2,il,ig)+ingrid(ir,ic-2,il,ig)-
     &    2.*ingrid(ir,ic,il,ig))/xincx**2
	 B21=B21*(-ingrid(ir+1,ic,il,it)+ingrid(ir-1,ic,il,it))/xincy
	 B21=B21*(ingrid(ir,ic,1,imap)**3.)
C Tal com estan aqui Q1 i Q2 necessitarem el signe de Fq al final, no l'hem
C inclos a dins Q en la seva expressio en termes del vent geostrofic
	  Q1(ir,ic,il)=g*(B20-B21)*1.e10/cor  

c CALCULO DE Q2
          B22=(ingrid(ir-2,ic,il,ig)+ingrid(ir+2,ic,il,ig)-
     &    2.*ingrid(ir,ic,il,ig))/xincy**2
          B22=B22*(ingrid(ir,ic+1,il,it)-ingrid(ir,ic-1,il,it))/xincx
	  B22=B22*(ingrid(ir,ic,1,imap)**3.)          
          B23=(-ingrid(ir+1,ic+1,il,ig)+ingrid(ir-1,ic+1,il,ig)+
     &	 ingrid(ir+1,ic-1,il,ig)-ingrid(ir-1,ic-1,il,ig))/(xincx*xincy)
	  B23=B23*(-ingrid(ir+1,ic,il,it)+ingrid(ir-1,ic,il,it))/xincy
	  B23=B23*(ingrid(ir,ic,1,imap)**3.)	  
	  Q2(ir,ic,il)=g*(B22-B23)*1.e10/cor
   
         enddo
        enddo 
C ASIGNACION DE LOS CONTORNOS DEL CAMPO DE Q1 Y Q2
        Do ir=3,nr-2
         Q1(ir,2,il)=2*Q1(ir,3,il)-Q1(ir,4,il)
         Q1(ir,1,il)=2*Q1(ir,2,il)-Q1(ir,3,il)
         Q1(ir,nc-1,il)=2*Q1(ir,nc-2,il)-Q1(ir,nc-3,il)
         Q1(ir,nc,il)=2*Q1(ir,nc-1,il)-Q1(ir,nc-2,il)
         Q2(ir,2,il)=2*Q2(ir,3,il)-Q2(ir,4,il)
         Q2(ir,1,il)=2*Q2(ir,2,il)-Q2(ir,3,il)
         Q2(ir,nc-1,il)=2*Q2(ir,nc-2,il)-Q2(ir,nc-3,il)
         Q2(ir,nc,il)=2*Q2(ir,nc-1,il)-Q2(ir,nc-2,il)
        Enddo

        Do ic=1,nc
         Q1(2,ic,il)=2*Q1(3,ic,il)-Q1(4,ic,il)
         Q1(1,ic,il)=2*Q1(2,ic,il)-Q1(3,ic,il)
         Q1(nr-1,ic,il)=2*Q1(nr-2,ic,il)-Q1(nr-3,ic,il)
         Q1(nr,ic,il)=2*Q1(nr-1,ic,il)-Q1(nr-2,ic,il)
         Q2(2,ic,il)=2*Q2(3,ic,il)-Q2(4,ic,il)
         Q2(1,ic,il)=2*Q2(2,ic,il)-Q2(3,ic,il)
         Q2(nr-1,ic,il)=2*Q2(nr-2,ic,il)-Q2(nr-3,ic,il)
         Q2(nr,ic,il)=2*Q2(nr-1,ic,il)-Q2(nr-2,ic,il)
        Enddo



c Calculo del forzamiento
C CALCULO DE LA DIVERGENCIA (DIVQ, en unidades de 1E-16 SI), ROTACIONAL 
C (ROTQ, en unidades de 1E-16 SI) Y FORZAMIENTO (FORZQ, en unidades de 
C 1E-18 SI)(Unitats de 1/(Pa s^3):


         do ic=2,nc-1
           do ir=2,nr-1
		lat=ingrid(ir,ic,1,ilat)*Pi/180.
		xincx=2.*ds
         DIV=(Q1(ir,ic+1,il)-Q1(ir,ic-1,il))/xincx+         
     $ (Q2(ir-1,ic,il)-Q2(ir+1,ic,il))/xincy
         DIV=DIV*ingrid(ir,ic,1,imap)     
C S´ha inclos el signe en aquesta expressio de Fq per que no l´haviem inclos
C a l´hora de calcular el camp Q
          FORZQ=-2.*(r/p(il))*DIV*1E8
	  outgrid(ir,ic,il) = forzQ
          enddo
        enddo 

C ASSIGNACION DE LOS CONTORNOS DEL CAMPO FORZQ

        do ir=2,nr-1
	 outgrid(ir,1,il)=2.*outgrid(ir,2,il)-outgrid(ir,3,il)
	 outgrid(ir,nc,il)=2.*outgrid(ir,nc-1,il)-outgrid(ir,nc-2,il)
	enddo
	do ic=2,nc-1
	outgrid(1,ic,il)=2.*outgrid(2,ic,il)-outgrid(3,ic,il)
	outgrid(nr,ic,il)=2.*outgrid(nr-1,ic,il)-outgrid(nr-2,ic,il)
	enddo
	outgrid(1,1,il)=2.*outgrid(2,2,il)-outgrid(3,3,il)
	outgrid(nr,nc,il)=2.*outgrid(nr-1,nc-1,il)-outgrid(nr-2,nc-2,il)
	outgrid(1,nc,il)=2.*outgrid(2,nc-1,il)-outgrid(3,nc-2,il)
	outgrid(nr,1,il)=2.*outgrid(nr-1,2,il)-outgrid(nr-2,3,il)

100      end do
      USERFUNC = 0
      RETURN

C  Bucle sobre todos los niveles de presion
143	print*,'Leyendo variables ya calculadas UQ y VQ'

        xincy=2.*ds
      do 105 il=1,maxnl
	p(il)=1012.5*exp(vert_args(il)/-7.2)*100.
c Calculo del forzamiento
C CALCULO DE LA DIVERGENCIA (DIVQ, en unidades de 1E-16 SI), ROTACIONAL 
C (ROTQ, en unidades de 1E-16 SI) Y FORZAMIENTO (FORZQ, en unidades de 
C 1E-18 SI):

         do ic=2,nc-1
           do ir=2,nr-1
		lat=(Proj_args(1)-(ir-1)*Proj_args(3))*Pi/180	
		xincx=2.*Rt*Cos(lat)*(Proj_args(4)*Pi/180.)
         DIV=(ingrid(ir,ic+1,il,iqu)-ingrid(ir,ic-1,il,iqu))/xincx        
     & +(ingrid(ir-1,ic,il,iqv)-ingrid(ir+1,ic,il,iqv))/xincy 
         DIV=DIV*ingrid(ir,ic,1,imap)           
          FORZQ=-2.*(r/p(il))*DIV*1E8
	  outgrid(ir,ic,il) = forzQ
          enddo
        enddo 
 
        do ir=2,nr-1
	 outgrid(ir,1,il)=2.*outgrid(ir,2,il)-outgrid(ir,3,il)
	 outgrid(ir,nc,il)=2.*outgrid(ir,nc-1,il)-outgrid(ir,nc-2,il)
	enddo
	do ic=2,nc-1
	outgrid(1,ic,il)=2.*outgrid(2,ic,il)-outgrid(3,ic,il)
	outgrid(nr,ic,il)=2.*outgrid(nr-1,ic,il)-outgrid(nr-2,ic,il)
	enddo
	outgrid(1,1,il)=2.*outgrid(2,2,il)-outgrid(3,3,il)
	outgrid(nr,nc,il)=2.*outgrid(nr-1,nc-1,il)-outgrid(nr-2,nc-2,il)
	outgrid(1,nc,il)=2.*outgrid(2,nc-1,il)-outgrid(3,nc-2,il)
	outgrid(nr,1,il)=2.*outgrid(nr-1,2,il)-outgrid(nr-2,3,il)

105      end do
      USERFUNC = 0
	print*,'Units of 10^-18 m/(kg s)'
      RETURN
      END
