*IF DEF,A15_1A                                                             CALCPV1A.2      
C ******************************COPYRIGHT******************************    GTS2F400.721    
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved.    GTS2F400.722    
C                                                                          GTS2F400.723    
C Use, duplication or disclosure of this code is subject to the            GTS2F400.724    
C restrictions as set forth in the contract.                               GTS2F400.725    
C                                                                          GTS2F400.726    
C                Meteorological Office                                     GTS2F400.727    
C                London Road                                               GTS2F400.728    
C                BRACKNELL                                                 GTS2F400.729    
C                Berkshire UK                                              GTS2F400.730    
C                RG12 2SZ                                                  GTS2F400.731    
C                                                                          GTS2F400.732    
C If no contract has been raised with this copy of the code, the use,      GTS2F400.733    
C duplication or disclosure of it is strictly prohibited.  Permission      GTS2F400.734    
C to do so must first be obtained in writing from the Head of Numerical    GTS2F400.735    
C Modelling at the above address.                                          GTS2F400.736    
C ******************************COPYRIGHT******************************    GTS2F400.737    
C                                                                          GTS2F400.738    
CLL Subroutine calc_pv -------------------------------------------------   CALCPV1A.3      
CLL                                                                        CALCPV1A.4      
CLL Purpose: To compute Ertel potential vorticity on isentropic            CALCPV1A.5      
CLL          surfaces.                                                     CALCPV1A.6      
CLL          Uses the Quasi-Hydrostatic equations, with complete           CALCPV1A.7      
CLL          representation of the Coriolis terms, and no metric           CALCPV1A.8      
CLL          terms omitted.                                                CALCPV1A.9      
CLL          The shallow atmosphere approximation is not made.             CALCPV1A.10     
CLL          Under UPDATE identifier GLOBAL, the data is                   CALCPV1A.11     
CLL          assumed periodic along the rows. Note that because            CALCPV1A.12     
CLL          it is a diagnostic routine care needs to be taken             CALCPV1A.13     
CLL          with missing data.                                            CALCPV1A.14     
CLL          Subroutine CALCRS1 is called to compute the psuedo radius.    CALCPV1A.15     
CLL          Values of p, rs, u, v and DTheta/Dp on theta surfaces         CALCPV1A.16     
CLL          are then computed assuming a cubic variation,                 CALCPV1A.17     
CLL          where possible.                                               CALCPV1A.18     
CLL          Two other subroutines then calculate the remaining            CALCPV1A.19     
CLL          terms in the potential vorticity expression.                  CALCPV1A.20     
CLL          The final part of the calc_pv routine then                    CALCPV1A.21     
CLL          adds in the rotational component and                          CALCPV1A.22     
CLL          checks for missing data.                                      CALCPV1A.23     
CLL                                                                        CALCPV1A.24     
CLL Not suitable for single column use.                                    CALCPV1A.25     
CLL                                                                        CALCPV1A.26     
CLL  7/10/92 Written by Simon Anderson                                     CALCPV1A.27     
CLL                                                                        CALCPV1A.28     
CLL  Model            Modification history from model version 3.0:         CALCPV1A.29     
CLL version  Date                                                          CALCPV1A.30     
CLL   3.1   14/01/93  Parameter list tidied up.                            MM180193.161    
CLL   3.2   28/07/93  Change subroutine name to uppercase for              TS280793.1      
CLL                   portability.    Tracey Smith                         TS280793.2      
CLL   3.4   26/05/94  Logical llints added to s/r args and passed to       GSS1F304.142    
CLL                                   calc_rs          S.J.Swarbrick       GSS1F304.143    
!LL   4.3   14/02/97  Added ARGFLDPT arguments and MPP code  P.Burton      GSM3F403.469    
CLL                                                                        CALCPV1A.31     
CLL Programming Standard: UM DOC Paper3, Version 4 (05/02/92)              CALCPV1A.32     
CLL                                                                        CALCPV1A.33     
CLL Logical Component Covered: D415                                        CALCPV1A.34     
CLL                                                                        CALCPV1A.35     
CLL System Task: D4                                                        CALCPV1A.36     
CLL                                                                        CALCPV1A.37     
CLL Documentation: U.M.D.P No 13. Derivation and Calculation of            CALCPV1A.38     
CLL                Unified Model Potential Vorticity.                      CALCPV1A.39     
CLL                by Simon Anderson and Ian Roulstone.                    CALCPV1A.40     
CLL                                                                        CALCPV1A.41     
CLLEND                                                                     CALCPV1A.42     
                                                                           CALCPV1A.43     
C*L ARGUMENTS: ---------------------------------------------------------   CALCPV1A.44     

      SUBROUTINE CALC_PV                                                    1,12TS280793.3      
     1                  (pstar,theta,u,v,p_field,u_field,                  CALCPV1A.46     
     2                   p_levels,row_length,                              CALCPV1A.47     
*CALL ARGFLDPT                                                             GSM3F403.470    
     3                   rmdi,ak,bk,des_theta,f3,                          CALCPV1A.48     
     & e_levels,n_levels,                                                  TD141293.58     
     & dth_dph,                                                            TD141293.59     
     4                   latitude_step_inverse,longitude_step_inverse,     CALCPV1A.49     
     5                   cos_u_latitude,sec_p_latitude,                    CALCPV1A.50     
     6                   pvort,llints)                                     GSS1F304.144    
                                                                           CALCPV1A.53     
                                                                           CALCPV1A.54     
      implicit none                                                        CALCPV1A.55     
      logical  llints                                                      GSS1F304.145    
                                                                           CALCPV1A.56     
C Input variables ------------------------------------------------------   CALCPV1A.57     
                                                                           CALCPV1A.58     
      integer                                                              CALCPV1A.59     
     & p_field                 !IN    Size of field on pressure points.    CALCPV1A.60     
     &,u_field                 !IN    Size of field on velocity points.    CALCPV1A.61     
     &,p_levels                !IN    Number of pressure levels.           CALCPV1A.62     
     &,row_length              !IN    Number of points in a row.           CALCPV1A.63     
     & ,n_levels       !IN Number of levels of dtheta/dp                   TD141293.60     
                                                                           CALCPV1A.64     
                                                                           GSM3F403.471    
*CALL TYPFLDPT                                                             GSM3F403.472    
      real                                                                 CALCPV1A.65     
     & pstar(p_field)          !IN    Surface pressure field.              CALCPV1A.66     
     &,u(u_field,p_levels)     !IN    Primary model array for u field.     CALCPV1A.67     
     &,v(u_field,p_levels)     !IN    Primary model array for v field.     CALCPV1A.68     
     &,theta(p_field,p_levels) !IN       "      "     "     theta field.   CALCPV1A.69     
                                                                           CALCPV1A.70     
      real                                                                 CALCPV1A.71     
     & rmdi                    !IN    Real missing data indicator.         CALCPV1A.72     
     &,ak(p_levels)            !IN    A coefficient of hybrid              CALCPV1A.73     
     &                         !      coordinates at full levels.          CALCPV1A.74     
     &,bk(p_levels)            !IN    B coefficient of hybrid              CALCPV1A.75     
     &                         !      coordinates at full levels.          CALCPV1A.76     
     & ,e_levels(n_levels)       !IN half-levels over range                TD141293.61     
     & ,dth_dph(p_field,n_levels)  !IN dtheta/dp half-levels               TD141293.62     
     &,des_theta               !IN    Value of theta we want pv on.        CALCPV1A.77     
     &,f3(u_field)             !IN    Coriolis term.                       CALCPV1A.78     
     &,latitude_step_inverse   !IN    1/latitude increment.                CALCPV1A.79     
     &,longitude_step_inverse  !IN    1/longitude increment.               CALCPV1A.80     
     &,cos_u_latitude(u_field) !IN    Cosine of latitude on u field.       CALCPV1A.81     
     &,sec_p_latitude(p_field) !IN    Secant of latitude on p field.       CALCPV1A.82     
                                                                           CALCPV1A.83     
                                                                           CALCPV1A.84     
C Output variables -----------------------------------------------------   CALCPV1A.85     
                                                                           CALCPV1A.86     
      real                                                                 CALCPV1A.87     
     & pvort(p_field)          !  OUT Value of potential vorticity         CALCPV1A.88     
     &                         !      on isentropic surface with           CALCPV1A.89     
     &                         !      theta=des_theta.                     CALCPV1A.90     
                                                                           CALCPV1A.91     
                                                                           CALCPV1A.97     
C*----------------------------------------------------------------------   CALCPV1A.98     
C*L Workspace Usage: 13 arrays are required.                               CALCPV1A.99     
      logical                                                              CALCPV1A.100    
     & mask(p_field)           ! Logical mask used to make                 CALCPV1A.101    
     &                         ! vectorising of a loop possible.           CALCPV1A.102    
      real                                                                 CALCPV1A.103    
     & rs(p_field,p_levels)    ! Calculated pseudo radius of Earth.        CALCPV1A.104    
     &,eta_level(p_field)    ! Interpolated eta-value of theta level       TD141293.63     
     &,u_on_theta(u_field)     ! Interpolated value of  u on theta surf.   CALCPV1A.106    
     &,u_p_on_theta(p_field)   ! Interpolated value of  u on theta surf.   CALCPV1A.107    
     &,v_on_theta(u_field)     ! Interpolated value of  v on theta surf.   CALCPV1A.108    
     &,v_p_on_theta(p_field)   ! Interpolated value of  v on theta surf.   CALCPV1A.109    
     &,rs_on_theta(p_field)    ! Interpolated value of rs on theta surf.   CALCPV1A.110    
     &,rs_uv_on_theta(u_field) ! Interpolated value of rs on theta surf.   CALCPV1A.111    
     &,dtheta_dp(p_field)      ! D(theta)/D(p) on theta surface.           CALCPV1A.112    
     &,vorticity1(p_field)     ! 1st Calculated term in the pv equation.   CALCPV1A.113    
     &,vorticity2(p_field)     ! 2nd Calculated term in the pv equation.   CALCPV1A.114    
     &,f3_p(p_field)           ! Interpolated f3 field on p field.         CALCPV1A.115    
                                                                           CALCPV1A.116    
C*----------------------------------------------------------------------   CALCPV1A.117    
C*L External subroutine calls:                                             CALCPV1A.118    
      external calc_rs         ! Compute pseudo radius.                    CALCPV1A.119    
      external pv_thint        ! Interpolate variables to theta level.     CALCPV1A.120    
      external vortic1         ! Compute the 1st vorticity term.           CALCPV1A.121    
      external vortic2         ! Compute the 2nd vorticity term.           CALCPV1A.122    
      external uv_to_p         ! Interpolate u grid field to p grid fld.   CALCPV1A.123    
                                                                           CALCPV1A.124    
C*----------------------------------------------------------------------   CALCPV1A.125    
C*L Call comdecks to get required variables:                               CALCPV1A.126    
*CALL C_G                                                                  CALCPV1A.127    
*CALL C_OMEGA                                                              CALCPV1A.128    
                                                                           CALCPV1A.129    
C*----------------------------------------------------------------------   CALCPV1A.130    
C*L Define local variables.                                                CALCPV1A.131    
      integer i,j,k          ! Loop variables.                             CALCPV1A.132    
      real mn                ! Mean value used in computing pole values.   CALCPV1A.133    
      real work1(p_field)    ! General workspace.                          CALCPV1A.134    
*IF DEF,MPP                                                                GSM3F403.473    
      INTEGER info !GCG return code                                        GSM3F403.474    
*IF -DEF,GLOBAL                                                            GSM3F403.475    
      INTEGER                                                              GSM3F403.476    
     &  START ! Work variables                                             GSM3F403.477    
     &  ,END                                                               GSM3F403.478    
*ENDIF                                                                     GSM3F403.479    
*ENDIF                                                                     GSM3F403.480    
                                                                           CALCPV1A.135    
                                                                           CALCPV1A.136    
C ----------------------------------------------------------------------   CALCPV1A.137    
CL Section 1 Compute rs, the pseudo radius.                                CALCPV1A.138    
CL ~~~~~~~~~                                                               CALCPV1A.139    
C ----------------------------------------------------------------------   CALCPV1A.140    
                                                                           CALCPV1A.141    
CL Section 1.1 Call CALC_RS to get rs for level 1.                         CALCPV1A.142    
CL ~~~~~~~~~~~ Rs is returned in rs(1,k).                                  CALCPV1A.143    
CL             Ts is returned in work1. Rs at level k-1 is input as        CALCPV1A.144    
CL             rs(1,2) as at k-1=0 the input is not used by calc_rs.       CALCPV1A.145    
                                                                           CALCPV1A.146    
      k=1                                                                  CALCPV1A.147    
      call calc_rs                                                         CALCPV1A.148    
     1            (pstar,ak,bk,work1,rs(1,2),rs(1,k),p_field,k,p_levels,   GSS1F304.146    
     2             llints)                                                 GSS1F304.147    
                                                                           CALCPV1A.150    
CL Section 1.2 Call CALC_RS to get rs for level k.                         CALCPV1A.151    
CL             Rs is returned in rs(1,k).                                  CALCPV1A.152    
CL             Ts is returned in work1. Rs at level k-1 is input as        CALCPV1A.153    
CL             rs(1,k-1).                                                  CALCPV1A.154    
CL             Loop from 2 to p_levels.                                    CALCPV1A.155    
                                                                           CALCPV1A.156    
      do 100 k=2,p_levels                                                  CALCPV1A.157    
                                                                           CALCPV1A.158    
        call calc_rs                                                       CALCPV1A.159    
     1          (pstar,ak,bk,work1,rs(1,k-1),rs(1,k),p_field,k,p_levels,   GSS1F304.148    
     2           llints)                                                   GSS1F304.149    
 100  continue                                                             CALCPV1A.162    
                                                                           CALCPV1A.163    
C ----------------------------------------------------------------------   CALCPV1A.164    
CL Section 2 Interpolate p, u, v and rs onto desired theta level,          CALCPV1A.165    
CL ~~~~~~~~~ and compute Dtheta/Dp assuming cubic variation if possible.   CALCPV1A.166    
C ----------------------------------------------------------------------   CALCPV1A.167    
                                                                           CALCPV1A.168    
      call pv_thint                                                        CALCPV1A.169    
     1             (pstar,theta,rs,u,v,p_field,u_field,                    CALCPV1A.170    
     2              p_levels,row_length,                                   CALCPV1A.171    
*CALL ARGFLDPT                                                             GSM3F403.481    
     3              rmdi,ak,bk,des_theta,                                  CALCPV1A.172    
     4 eta_level,rs_on_theta,rs_uv_on_theta,                               TD141293.64     
     & e_levels,dth_dph,n_levels,                                          TD141293.65     
     5              u_on_theta,v_on_theta,dtheta_dp)                       CALCPV1A.174    
                                                                           CALCPV1A.175    
                                                                           CALCPV1A.176    
C ----------------------------------------------------------------------   CALCPV1A.177    
CL Section 3 Compute the various terms in the potential vorticity          CALCPV1A.178    
CL ~~~~~~~~~ equation for desired theta surface.                           CALCPV1A.179    
C ----------------------------------------------------------------------   CALCPV1A.180    
                                                                           CALCPV1A.181    
CL Section 3.1 Compute -2*omega*cos(phi)/rs D(rs)/D(phi).                  CALCPV1A.182    
CL ~~~~~~~~~~~                                                             CALCPV1A.183    
      call vortic1                                                         CALCPV1A.184    
     1            (rs_on_theta,                                            CALCPV1A.185    
     2             sec_p_latitude,vorticity1,                              CALCPV1A.186    
     3             p_field,row_length,                                     CALCPV1A.187    
*CALL ARGFLDPT                                                             GSM3F403.482    
     4             latitude_step_inverse,longitude_step_inverse)           CALCPV1A.188    
                                                                           CALCPV1A.189    
                                                                           CALCPV1A.190    
CL Section 3.2 Compute 1/(rs*rs*cos(phi)) *                                CALCPV1A.191    
CL ~~~~~~~~~~~         (D(rs.v)/D(lambda)-D(rs.cos(phi).u)/D(phi)).        CALCPV1A.192    
      call vortic2                                                         CALCPV1A.193    
     1            (u_on_theta,v_on_theta,rs_on_theta,                      CALCPV1A.194    
     2             rs_uv_on_theta,cos_u_latitude,                          CALCPV1A.195    
     3             sec_p_latitude,vorticity2,                              CALCPV1A.196    
     4             p_field,u_field,row_length,                             CALCPV1A.197    
*CALL ARGFLDPT                                                             GSM3F403.483    
     5             latitude_step_inverse,longitude_step_inverse)           CALCPV1A.198    
                                                                           CALCPV1A.199    
                                                                           CALCPV1A.200    
C ----------------------------------------------------------------------   CALCPV1A.201    
CL Section 4 Compute the potential vorticity using the vorticity terms,    CALCPV1A.202    
CL ~~~~~~~~~ D(theta)/D(p) and the value of the coriolis term.             CALCPV1A.203    
CL           Care needs to be taken for rmdi.                              CALCPV1A.204    
C ----------------------------------------------------------------------   CALCPV1A.205    
                                                                           CALCPV1A.206    
CL Section 4.1 Firstly, find which points we can actually calculate        CALCPV1A.207    
CL ~~~~~~~~~~~ potential vorticity on. To do this we test for missing      CALCPV1A.208    
CL             data at all points in a 3x3 stencil around each point.      CALCPV1A.209    
CL             Special care is taken at the boundaries and poles.          CALCPV1A.210    
                                                                           CALCPV1A.211    
C Set logical mask default to .TRUE.                                       CALCPV1A.212    
      do i = FIRST_FLD_PT,LAST_P_FLD_PT                                    GSM3F403.484    
        mask(i) = .true.                                                   CALCPV1A.214    
      end do                                                               CALCPV1A.215    
                                                                           CALCPV1A.216    
C Check points on non-polar rows:                                          CALCPV1A.217    
      do j=FIRST_ROW,P_LAST_ROW                                            GSM3F403.485    
C First point on each row:                                                 CALCPV1A.219    
        i=(j-1)*row_length+FIRST_ROW_PT                                    GSM3F403.486    
*IF DEF,GLOBAL                                                             CALCPV1A.221    
*IF -DEF,MPP                                                               GSM3F403.487    
        if (   eta_level(i-1           ).eq.rmdi.                          TD141293.66     
     &      or.eta_level(i-row_length  ).eq.rmdi.                          TD141293.67     
     &      or.eta_level(i-row_length+1).eq.rmdi.                          TD141293.68     
     &      or.eta_level(i+row_length-1).eq.rmdi.                          TD141293.69     
     &      or.eta_level(i             ).eq.rmdi.                          TD141293.70     
     &      or.eta_level(i+1           ).eq.rmdi.                          TD141293.71     
     &      or.eta_level(i+2*row_length-1).eq.rmdi.                        TD141293.72     
     &      or.eta_level(i+row_length  ).eq.rmdi.                          TD141293.73     
     &      or.eta_level(i+row_length+1).eq.rmdi)   mask(i)=.false.        TD141293.74     
*ENDIF                                                                     GSM3F403.488    
*ELSE                                                                      GSM3F403.489    
*IF DEF,MPP                                                                GSM3F403.490    
        IF (at_left_of_LPG) THEN                                           GSM3F403.491    
          mask(i+EW_Halo)=.false.                                          GSM3F403.492    
        ENDIF                                                              GSM3F403.493    
*ELSE                                                                      CALCPV1A.231    
        mask(i)=.false.                                                    CALCPV1A.232    
*ENDIF                                                                     CALCPV1A.233    
*ENDIF                                                                     GSM3F403.494    
                                                                           GSM3F403.495    
C Inner points on each row:                                                CALCPV1A.234    
*IF DEF,MPP                                                                GSM3F403.496    
*IF -DEF,GLOBAL                                                            GSM3F403.497    
C Borders of area done separately                                          GSM3F403.498    
        START=EW_Halo+1                                                    GSM3F403.499    
        END=EW_Halo                                                        GSM3F403.500    
        IF (at_left_of_LPG) START=START+1                                  GSM3F403.501    
        IF (at_right_of_LPG) END=END+1                                     GSM3F403.502    
        do i=(j-1)*row_length+START,(j-1)*row_length+row_length-END        GSM3F403.503    
*ELSE                                                                      GSM3F403.504    
        do i=(j-1)*row_length+1+EW_Halo,                                   GSM3F403.505    
     &      j*row_length-EW_Halo                                           GSM3F403.506    
*ENDIF                                                                     GSM3F403.507    
*ELSE                                                                      GSM3F403.508    
        do i=(j-1)*row_length+2,(j-1)*row_length+row_length-1              CALCPV1A.235    
*ENDIF                                                                     GSM3F403.509    
        if (   eta_level(i-row_length-1).eq.rmdi.                          TD141293.75     
     &      or.eta_level(i-row_length  ).eq.rmdi.                          TD141293.76     
     &      or.eta_level(i-row_length+1).eq.rmdi.                          TD141293.77     
     &      or.eta_level(i-1           ).eq.rmdi.                          TD141293.78     
     &      or.eta_level(i             ).eq.rmdi.                          TD141293.79     
     &      or.eta_level(i+1           ).eq.rmdi.                          TD141293.80     
     &      or.eta_level(i+row_length-1).eq.rmdi.                          TD141293.81     
     &      or.eta_level(i+row_length  ).eq.rmdi.                          TD141293.82     
     &      or.eta_level(i+row_length+1).eq.rmdi)   mask(i)=.false.        TD141293.83     
        end do                                                             CALCPV1A.245    
C Last point on each row:                                                  CALCPV1A.246    
        i=(j-1)*row_length+LAST_ROW_PT                                     GSM3F403.510    
*IF DEF,GLOBAL                                                             CALCPV1A.248    
*IF -DEF,MPP                                                               GSM3F403.511    
        if (   eta_level(i-row_length-1).eq.rmdi.                          TD141293.84     
     &      or.eta_level(i-row_length  ).eq.rmdi.                          TD141293.85     
     &      or.eta_level(i-2*row_length+1).eq.rmdi.                        TD141293.86     
     &      or.eta_level(i-1           ).eq.rmdi.                          TD141293.87     
     &      or.eta_level(i             ).eq.rmdi.                          TD141293.88     
     &      or.eta_level(i-row_length+1).eq.rmdi.                          TD141293.89     
     &      or.eta_level(i+row_length-1).eq.rmdi.                          TD141293.90     
     &      or.eta_level(i+row_length  ).eq.rmdi.                          TD141293.91     
     &      or.eta_level(i+1           ).eq.rmdi)   mask(i)=.false.        TD141293.92     
*ENDIF                                                                     GSM3F403.512    
*ELSE                                                                      GSM3F403.513    
*IF DEF,MPP                                                                GSM3F403.514    
        IF (at_right_of_LPG) THEN                                          GSM3F403.515    
          mask(i-EW_Halo)=.false.                                          GSM3F403.516    
        ENDIF                                                              GSM3F403.517    
*ELSE                                                                      CALCPV1A.258    
        mask(i)=.false.                                                    CALCPV1A.259    
*ENDIF                                                                     GSM3F403.518    
*ENDIF                                                                     CALCPV1A.260    
      end do                                                               CALCPV1A.261    
                                                                           CALCPV1A.262    
*IF DEF,GLOBAL                                                             CALCPV1A.263    
C Check North pole:                                                        CALCPV1A.264    
      j=0                                                                  CALCPV1A.265    
*IF DEF,MPP                                                                GSM3F403.519    
      IF (at_top_of_LPG) THEN                                              GSM3F403.520    
*ENDIF                                                                     GSM3F403.521    
      if (eta_level(TOP_ROW_START+LAST_ROW_PT-1) .eq. rmdi) then           GSM3F403.522    
        j=1                                                                GSM3F403.523    
      ELSE                                                                 GSM3F403.524    
        do i=TOP_ROW_START+ROW_LENGTH,                                     GSM3F403.525    
     &       TOP_ROW_START+ROW_LENGTH+LAST_ROW_PT-1                        GSM3F403.526    
          if(eta_level(i).eq.rmdi) j=1                                     GSM3F403.527    
        enddo                                                              GSM3F403.528    
      endif                                                                GSM3F403.529    
*IF DEF,MPP                                                                GSM3F403.530    
      CALL GCG_IMAX(1,GC_ROW_GROUP,info,j)                                 GSM3F403.531    
*ENDIF                                                                     GSM3F403.532    
      if (j.eq.1) then                                                     GSM3F403.533    
        do i=TOP_ROW_START,TOP_ROW_START+LAST_ROW_PT-1                     GSM3F403.534    
          mask(i)=.false.                                                  GSM3F403.535    
        enddo                                                              GSM3F403.536    
      endif                                                                GSM3F403.537    
*IF DEF,MPP                                                                GSM3F403.538    
      endif ! if at_top_of_LPG                                             GSM3F403.539    
*ENDIF                                                                     GSM3F403.540    
                                                                           CALCPV1A.273    
C Check South pole:                                                        CALCPV1A.274    
      j=0                                                                  CALCPV1A.275    
*IF DEF,MPP                                                                GSM3F403.541    
      IF (at_base_of_LPG) THEN                                             GSM3F403.542    
*ENDIF                                                                     GSM3F403.543    
      if (eta_level(P_BOT_ROW_START) .eq. rmdi) then                       GSM3F403.544    
        j=1                                                                GSM3F403.545    
      else                                                                 GSM3F403.546    
        do i=P_BOT_ROW_START-ROW_LENGTH,                                   GSM3F403.547    
     &       P_BOT_ROW_START-ROW_LENGTH+LAST_ROW_PT-1                      GSM3F403.548    
          if (eta_level(i) .eq. rmdi) j=1                                  GSM3F403.549    
        enddo                                                              GSM3F403.550    
      endif                                                                GSM3F403.551    
*IF DEF,MPP                                                                GSM3F403.552    
      CALL GCG_IMAX(1,GC_ROW_GROUP,info,j)                                 GSM3F403.553    
*ENDIF                                                                     GSM3F403.554    
                                                                           GSM3F403.555    
      if (j.eq.1) then                                                     GSM3F403.556    
        do i=P_BOT_ROW_START,P_BOT_ROW_START+LAST_ROW_PT-1                 GSM3F403.557    
          mask(i)=.false.                                                  GSM3F403.558    
        enddo                                                              GSM3F403.559    
      endif                                                                GSM3F403.560    
*IF DEF,MPP                                                                GSM3F403.561    
      endif ! if at_base_of_LPG                                            GSM3F403.562    
*ENDIF                                                                     GSM3F403.563    
*ELSE                                                                      CALCPV1A.283    
C Set mask on Northern and Southern boundaries:                            CALCPV1A.284    
*IF DEF,MPP                                                                GSM3F403.564    
      IF (at_top_of_LPG) THEN                                              GSM3F403.565    
*ENDIF                                                                     GSM3F403.566    
        do i=TOP_ROW_START,TOP_ROW_START+LAST_ROW_PT-1                     GSM3F403.567    
          mask(i)=.false.                                                  GSM3F403.568    
        enddo                                                              GSM3F403.569    
*IF DEF,MPP                                                                GSM3F403.570    
      ENDIF ! if at_top_of_LPG                                             GSM3F403.571    
                                                                           GSM3F403.572    
      IF (at_base_of_LPG) THEN                                             GSM3F403.573    
*ENDIF                                                                     GSM3F403.574    
        do i=P_BOT_ROW_START,P_BOT_ROW_START+LAST_ROW_PT-1                 GSM3F403.575    
          mask(i)=.false.                                                  GSM3F403.576    
        enddo                                                              GSM3F403.577    
*IF DEF,MPP                                                                GSM3F403.578    
      ENDIF ! if at_base_of_LPG                                            GSM3F403.579    
*ENDIF                                                                     GSM3F403.580    
*ENDIF                                                                     GSM3F403.581    
*IF DEF,MPP                                                                GSM3F403.582    
! Set all the halos to false                                               GSM3F403.583    
      if (at_top_of_lpg) then                                              GSM3F403.584    
        do j=1,NS_Halo                                                     GSM3F403.585    
          do i=(j-1)*row_length+1,j*row_length                             GSM3F403.586    
            mask(i)=.false.                                                GSM3F403.587    
          enddo                                                            GSM3F403.588    
        enddo                                                              GSM3F403.589    
      endif ! if at_top_of_lpg                                             GSM3F403.590    
                                                                           GSM3F403.591    
      if (at_base_of_lpg) then                                             GSM3F403.592    
        do j=tot_P_ROWS-NS_Halo+1,tot_P_ROWS                               GSM3F403.593    
          do i=(j-1)*row_length+1,j*row_length                             GSM3F403.594    
            mask(i)=.false.                                                GSM3F403.595    
          enddo                                                            GSM3F403.596    
        enddo                                                              GSM3F403.597    
      endif ! if at_base_of_lpg                                            GSM3F403.598    
                                                                           GSM3F403.599    
      if (at_left_of_lpg) THEN                                             GSM3F403.600    
        do j=NS_Halo+1,tot_P_ROWS-NS_Halo                                  GSM3F403.601    
          do i=(j-1)*row_length+1,(j-1)*row_length+EW_Halo                 GSM3F403.602    
            mask(i)=.false.                                                GSM3F403.603    
          enddo                                                            GSM3F403.604    
        enddo                                                              GSM3F403.605    
      endif ! if at_left_of_lpg                                            GSM3F403.606    
                                                                           GSM3F403.607    
      if (at_right_of_lpg) THEN                                            GSM3F403.608    
        do j=NS_Halo+1,tot_P_ROWS-NS_Halo                                  GSM3F403.609    
          do i=j*row_length-EW_Halo+1,j*row_length                         GSM3F403.610    
            mask(i)=.false.                                                GSM3F403.611    
          enddo                                                            GSM3F403.612    
        enddo                                                              GSM3F403.613    
      endif                                                                GSM3F403.614    
*ENDIF                                                                     CALCPV1A.291    
                                                                           CALCPV1A.292    
                                                                           CALCPV1A.293    
CL Section 4.2 Interpolate fields u, v, and F3 (the Coriolis term)         CALCPV1A.294    
CL ~~~~~~~~~~~ to the p-grid.                                              CALCPV1A.295    
                                                                           CALCPV1A.296    
*IF DEF,MPP                                                                GSM3F403.615    
      call uv_to_p                                                         GSM3F403.616    
     1           (u_on_theta,u_p_on_theta(row_length+1),u_field,           GSM3F403.617    
     2            p_field-row_length,row_length,U_LAST_ROW+1)              GSM3F403.618    
      call uv_to_p                                                         GSM3F403.619    
     1           (v_on_theta,v_p_on_theta(row_length+1),u_field,           GSM3F403.620    
     2            p_field-row_length,row_length,U_LAST_ROW+1)              GSM3F403.621    
      call uv_to_p                                                         GSM3F403.622    
     1           (f3,f3_p(row_length+1),u_field,                           GSM3F403.623    
     2            p_field-row_length,row_length,U_LAST_ROW+1)              GSM3F403.624    
*ELSE                                                                      GSM3F403.625    
      call uv_to_p                                                         CALCPV1A.297    
     1           (u_on_theta,u_p_on_theta(row_length+1),u_field,p_field,   CALCPV1A.298    
     2            row_length,u_field/row_length)                           CALCPV1A.299    
      call uv_to_p                                                         CALCPV1A.300    
     1           (v_on_theta,v_p_on_theta(row_length+1),u_field,p_field,   CALCPV1A.301    
     2            row_length,u_field/row_length)                           CALCPV1A.302    
      call uv_to_p                                                         CALCPV1A.303    
     1           (f3,f3_p(row_length+1),u_field,p_field,                   CALCPV1A.304    
     2            row_length,u_field/row_length)                           CALCPV1A.305    
*ENDIF                                                                     GSM3F403.626    
                                                                           CALCPV1A.307    
CL Section 4.3 Calculate the potential vorticity.                          CALCPV1A.308    
CL ~~~~~~~~~~~                                                             CALCPV1A.309    
                                                                           CALCPV1A.310    
      do i = START_POINT_NO_HALO,END_P_POINT_NO_HALO                       GSM3F403.627    
        if (mask(i)) then                                                  CALCPV1A.312    
                                                                           CALCPV1A.313    
          pvort(i) = (vorticity1(i) + f3_p(i) + vorticity2(i)) *           CALCPV1A.314    
     &               dtheta_dp(i) *                                        CALCPV1A.315    
     &               ((2.*omega + u_p_on_theta(i)*                         CALCPV1A.316    
     &                 sec_p_latitude(i)/rs_on_theta(i))*                  CALCPV1A.317    
     &               (u_p_on_theta(i)/sec_p_latitude(i)) +                 CALCPV1A.318    
     &               (v_p_on_theta(i)*v_p_on_theta(i)/rs_on_theta(i))-g)   CALCPV1A.319    
                                                                           CALCPV1A.320    
        else                                                               CALCPV1A.321    
          pvort(i) = rmdi                                                  CALCPV1A.322    
        end if                                                             CALCPV1A.323    
      enddo                                                                CALCPV1A.324    
                                                                           CALCPV1A.325    
CL Section 4.3.1 Compute the potential vorticity at the Northern           CALCPV1A.326    
CL ~~~~~~~~~~~~~ and Southern boundaries.                                  CALCPV1A.327    
CL               Note that under UPDATE identifier GLOBAL, as pv           CALCPV1A.328    
CL               is a scaler this is defined to be the mean of             CALCPV1A.329    
CL               the near pole row.                                        CALCPV1A.330    
CL               If any rmdi is present in the row next to the pole,       CALCPV1A.331    
CL               then the pole has a value rmdi.                           CALCPV1A.332    
CL               ELSE, the value returned is missing data.                 CALCPV1A.333    
                                                                           CALCPV1A.334    
C Work out the Northern boundary.                                          CALCPV1A.335    
*IF DEF,MPP                                                                GSM3F403.628    
      if (at_top_of_LPG) then                                              GSM3F403.629    
*ENDIF                                                                     GSM3F403.630    
        if (mask(TOP_ROW_START+LAST_ROW_PT-1)) then                        GSM3F403.631    
                                                                           CALCPV1A.352    
          mn=0                                                             GSM3F403.632    
*IF -DEF,MPP                                                               GSM3F403.633    
          do i=TOP_ROW_START+row_length,                                   GSM3F403.634    
     &         TOP_ROW_START+2*row_length-1                                GSM3F403.635    
            mn=mn+pvort(i)                                                 GSM3F403.636    
          enddo                                                            GSM3F403.637    
*ELSE                                                                      GSM3F403.638    
*IF DEF,REPROD                                                             GSM3F403.639    
          call GCG_RVECSUMR(                                               GSM3F403.640    
*ELSE                                                                      GSM3F403.641    
          call GCG_RVECSUMF(                                               GSM3F403.642    
*ENDIF                                                                     GSM3F403.643    
     &      row_length-2*EW_Halo,row_length-2*EW_Halo,1,1,                 GSM3F403.644    
     &      pvort(TOP_ROW_START+FIRST_ROW_PT-1+row_length),                GSM3F403.645    
     &      GC_ROW_GROUP,info,mn)                                          GSM3F403.646    
*ENDIF                                                                     GSM3F403.647    
          mn=mn/GLOBAL_ROW_LENGTH                                          GSM3F403.648    
                                                                           GSM3F403.649    
          do i=TOP_ROW_START,TOP_ROW_START+row_length-1                    GSM3F403.650    
            pvort(i)=mn                                                    GSM3F403.651    
          enddo                                                            GSM3F403.652    
        else                                                               GSM3F403.653    
! Missing data                                                             GSM3F403.654    
          do i=TOP_ROW_START,TOP_ROW_START+row_length-1                    GSM3F403.655    
            pvort(i)=rmdi                                                  GSM3F403.656    
          enddo                                                            GSM3F403.657    
        endif                                                              GSM3F403.658    
*IF DEF,MPP                                                                GSM3F403.659    
      endif ! if at_top_of_LPG                                             GSM3F403.660    
*ENDIF                                                                     GSM3F403.661    
                                                                           GSM3F403.662    
! Southern boundary                                                        GSM3F403.663    
*IF DEF,MPP                                                                GSM3F403.664    
                                                                           GSM3F403.665    
      if (at_base_of_LPG) then                                             GSM3F403.666    
*ENDIF                                                                     GSM3F403.667    
        if (mask(P_BOT_ROW_START+FIRST_ROW_PT-1)) then                     GSM3F403.668    
                                                                           GSM3F403.669    
          mn=0                                                             GSM3F403.670    
*IF -DEF,MPP                                                               GSM3F403.671    
          do i=P_BOT_ROW_START-ROW_LENGTH,                                 GSM3F403.672    
     &         P_BOT_ROW_START-1                                           GSM3F403.673    
            mn=mn+pvort(i)                                                 GSM3F403.674    
          enddo                                                            GSM3F403.675    
*ELSE                                                                      GSM3F403.676    
*IF DEF,REPROD                                                             GSM3F403.677    
          call GCG_RVECSUMR(                                               GSM3F403.678    
*ELSE                                                                      GSM3F403.679    
          call GCG_RVECSUMF(                                               GSM3F403.680    
*ENDIF                                                                     GSM3F403.681    
     &      row_length-2*EW_Halo,row_length-2*EW_Halo,1,1,                 GSM3F403.682    
     &      pvort(P_BOT_ROW_START+FIRST_ROW_PT-1-row_length),              GSM3F403.683    
     &      GC_ROW_GROUP,info,mn)                                          GSM3F403.684    
*ENDIF                                                                     GSM3F403.685    
          mn=mn/GLOBAL_ROW_LENGTH                                          GSM3F403.686    
                                                                           GSM3F403.687    
          do i=P_BOT_ROW_START,P_BOT_ROW_START+row_length-1                GSM3F403.688    
            pvort(i)=mn                                                    GSM3F403.689    
          enddo                                                            GSM3F403.690    
        else                                                               GSM3F403.691    
! Missing data                                                             GSM3F403.692    
          do i=P_BOT_ROW_START,P_BOT_ROW_START+row_length-1                GSM3F403.693    
            pvort(i)=rmdi                                                  GSM3F403.694    
          enddo                                                            GSM3F403.695    
        endif                                                              GSM3F403.696    
*IF DEF,MPP                                                                GSM3F403.697    
      endif ! if at_base_of_LPG                                            GSM3F403.698    
*ENDIF                                                                     GSM3F403.699    
                                                                           GSM3F403.700    
*IF DEF,MPP                                                                GSM3F403.701    
! Fill in the N+S halo areas with rmdi                                     GSM3F403.702    
      do i=1,FIRST_FLD_PT-1                                                GSM3F403.703    
        pvort(i)=rmdi                                                      GSM3F403.704    
      enddo                                                                GSM3F403.705    
      do i=LAST_P_FLD_PT+1,p_field                                         GSM3F403.706    
        pvort(i)=rmdi                                                      GSM3F403.707    
      enddo                                                                GSM3F403.708    
      CALL SWAPBOUNDS(pvort,ROW_LENGTH,tot_P_ROWS,EW_Halo,                 GSM3F403.709    
     &  NS_Halo,1)                                                         GSM3F403.710    
*ENDIF                                                                     GSM3F403.711    
                                                                           CALCPV1A.370    
                                                                           CALCPV1A.371    
      return                                                               CALCPV1A.372    
      end                                                                  CALCPV1A.373    
                                                                           CALCPV1A.374    
*ENDIF                                                                     CALCPV1A.375