*IF DEF,A15_1A                                                             CALCP21A.2      
C ******************************COPYRIGHT******************************    GTS2F400.703    
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved.    GTS2F400.704    
C                                                                          GTS2F400.705    
C Use, duplication or disclosure of this code is subject to the            GTS2F400.706    
C restrictions as set forth in the contract.                               GTS2F400.707    
C                                                                          GTS2F400.708    
C                Meteorological Office                                     GTS2F400.709    
C                London Road                                               GTS2F400.710    
C                BRACKNELL                                                 GTS2F400.711    
C                Berkshire UK                                              GTS2F400.712    
C                RG12 2SZ                                                  GTS2F400.713    
C                                                                          GTS2F400.714    
C If no contract has been raised with this copy of the code, the use,      GTS2F400.715    
C duplication or disclosure of it is strictly prohibited.  Permission      GTS2F400.716    
C to do so must first be obtained in writing from the Head of Numerical    GTS2F400.717    
C Modelling at the above address.                                          GTS2F400.718    
C ******************************COPYRIGHT******************************    GTS2F400.719    
C                                                                          GTS2F400.720    
CLL Subroutine calc_pv_p -----------------------------------------------   CALCP21A.3      
CLL                                                                        CALCP21A.4      
CLL Purpose: To compute Ertel potential vorticity                          CALCP21A.5      
CLL          on pressure levels.                                           CALCP21A.6      
CLL          Uses the Quasi-Hydrostatic equations, with complete           CALCP21A.7      
CLL          representation of the Coriolis terms, and no metric           CALCP21A.8      
CLL          terms omitted.                                                CALCP21A.9      
CLL          The shallow atmosphere approximation is not made.             CALCP21A.10     
CLL          Under UPDATE identifier GLOBAL, the data is                   CALCP21A.11     
CLL          assumed periodic along the rows. Note that because            CALCP21A.12     
CLL          it is a diagnostic routine care needs to be taken             CALCP21A.13     
CLL          with missing data.                                            CALCP21A.14     
CLL                                                                        CALCP21A.15     
CLL Not suitable for single column use.                                    CALCP21A.16     
CLL                                                                        CALCP21A.17     
CLL  Model            Modification history:                                CALCP21A.18     
CLL Version   Date                                                         CALCP21A.19     
CLL   3.1   12/11/92  Written by Simon Anderson.                           CALCP21A.20     
CLL   3.1   18/01/93  New deck at the release of Version 3.1.              CALCP21A.21     
CLL   3.2   28/07/93  Change subroutine name to uppercase for              TS280793.4      
CLL                   portability.    Tracey Smith                         TS280793.5      
CLL   3.4   26/05/94  Argument llints added and passed to calc_rs          GSS1F304.150    
CLL                                                     S.J.Swarbrick      GSS1F304.151    
CLL                                                                        CALCP21A.22     
CLL Programming Standard: UM DOC Paper3, Version 4 (05/02/92)              CALCP21A.23     
CLL                                                                        CALCP21A.24     
CLL Logical Component Covered: D415                                        CALCP21A.25     
CLL                                                                        CALCP21A.26     
CLL System Task: D4                                                        CALCP21A.27     
CLL                                                                        CALCP21A.28     
CLL Documentation: U.M.D.P No.13. Derivation and Calculation of            CALCP21A.29     
CLL                Unified Model Potential Vorticity.                      CALCP21A.30     
CLL                By Simon Anderson and Ian Roulstone.                    CALCP21A.31     
CLL                                                                        CALCP21A.32     
CLLEND                                                                     CALCP21A.33     
C                                                                          CALCP21A.34     
C*L ARGUMENTS: ---------------------------------------------------------   CALCP21A.35     

      SUBROUTINE CALC_PV_P                                                  2,9TS280793.6      
     1                    (pstar,theta,u,v,p_field,u_field,                CALCP21A.37     
     2                     p_levels,row_length,                            CALCP21A.38     
*CALL ARGFLDPT                                                             GSM3F403.712    
     3                     rmdi,ak,bk,des_press,f3,                        CALCP21A.39     
     & e_levels,n_levels,dth_dph,                                          TD141293.175    
     4                     latitude_step_inverse,longitude_step_inverse,   CALCP21A.40     
     5                     cos_u_latitude,sec_p_latitude,                  CALCP21A.41     
     6                     pvort_p,theta_on_press,llints)                  GSS1F304.152    
                                                                           CALCP21A.43     
                                                                           CALCP21A.44     
      implicit none                                                        CALCP21A.45     
      logical  llints                                                      GSS1F304.153    
                                                                           CALCP21A.46     
C Input variables ------------------------------------------------------   CALCP21A.47     
                                                                           CALCP21A.48     
      integer                                                              CALCP21A.49     
     & p_field                 !IN    Size of field on pressure points.    CALCP21A.50     
     &,u_field                 !IN    Size of field on velocity points.    CALCP21A.51     
     &,p_levels                !IN    Number of pressure levels.           CALCP21A.52     
     &,row_length              !IN    Number of points in a row.           CALCP21A.53     
     & ,n_levels           !IN Number of levels of spline                  TD141293.176    
                                                                           CALCP21A.54     
*CALL TYPFLDPT                                                             GSM3F403.713    
                                                                           GSM3F403.714    
      real                                                                 CALCP21A.55     
     & pstar(p_field)          !IN    Surface pressure field.              CALCP21A.56     
     &,u(u_field,p_levels)     !IN    Primary model array for u field.     CALCP21A.57     
     &,v(u_field,p_levels)     !IN    Primary model array for v field.     CALCP21A.58     
     &,theta(p_field,p_levels) !IN       "      "     "     theta field.   CALCP21A.59     
     & ,dth_dph(p_field,n_levels)  !IN dtheta/dp half-levels               TD141293.177    
                                                                           CALCP21A.60     
      real                                                                 CALCP21A.61     
     & rmdi                    !IN    Real missing data indicator.         CALCP21A.62     
     &,ak(p_levels)            !IN    A coefficient of hybrid              CALCP21A.63     
     &                         !      coordinates at full levels.          CALCP21A.64     
     &,bk(p_levels)            !IN    B coefficient of hybrid              CALCP21A.65     
     &                         !      coordinates at full levels.          CALCP21A.66     
     & ,e_levels(n_levels)       !IN half-levels over range                TD141293.178    
     &,des_press               !IN    Value of pressure we want pv on.     CALCP21A.67     
     &,f3(u_field)             !IN    Coriolis term.                       CALCP21A.68     
     &,latitude_step_inverse   !IN    1/latitude increment.                CALCP21A.69     
     &,longitude_step_inverse  !IN    1/longitude increment.               CALCP21A.70     
     &,cos_u_latitude(u_field) !IN    Cosine of latitude on u field.       CALCP21A.71     
     &,sec_p_latitude(p_field) !IN    Secant of latitude on p field.       CALCP21A.72     
                                                                           CALCP21A.73     
                                                                           CALCP21A.74     
C Output variables -----------------------------------------------------   CALCP21A.75     
                                                                           CALCP21A.76     
      real                                                                 CALCP21A.77     
     & pvort_p(p_field)        !  OUT Value of potential vorticity         CALCP21A.78     
     &                         !      on pressure level with               CALCP21A.79     
     &                         !      pressure=des_press.                  CALCP21A.80     
     &,theta_on_press(p_field) !  OUT Value of theta on pressure           CALCP21A.81     
     &                         !      level with pressure=des_press.       CALCP21A.82     
                                                                           CALCP21A.83     
C*----------------------------------------------------------------------   CALCP21A.84     
C*L Workspace Usage: 14 arrays are required.                               CALCP21A.85     
      logical                                                              CALCP21A.86     
     & mask(p_field)           ! Logical mask used to make                 CALCP21A.87     
     &                         ! vectorising of a loop possible.           CALCP21A.88     
      real                                                                 CALCP21A.89     
     & rs(p_field,p_levels)    ! Calculated pseudo radius of Earth.        CALCP21A.90     
     &,u_on_press(u_field)     ! Interpolated value of  u on p level.      CALCP21A.91     
     &,v_on_press(u_field)     ! Interpolated value of  v on p level.      CALCP21A.92     
     &,rs_on_press(p_field)    ! Interpolated value of rs on p level.      CALCP21A.93     
     &,drsu_dp(p_field)        ! D(rs.u)/D(p) on pressure level.           CALCP21A.94     
     &,drsv_dp(p_field)        ! D(rs.v)/D(p) on pressure level.           CALCP21A.95     
     &,dtheta_dp(p_field)      ! D(Theta)/D(p) on pressure level.          CALCP21A.96     
     &,dtheta_dlatitude(p_field) ! D(Theta)/D(lat) on pressure level.      CALCP21A.97     
     &,vorticity3(p_field)     ! A calculated term in the pv equation.     CALCP21A.98     
     &,vorticity4(p_field)     ! A calculated term in the pv equation.     CALCP21A.99     
     &,vorticity5(p_field)     ! A calculated term in the pv equation.     CALCP21A.100    
     &,f3_p(p_field)           ! Interpolated f3 field on p field.         CALCP21A.101    
     &,eta_level(p_field)    ! Interpolated eta-value of theta level       TD141293.179    
     &,work1(p_field)          ! General workspace.                        CALCP21A.102    
                                                                           CALCP21A.103    
C*----------------------------------------------------------------------   CALCP21A.104    
C*L External subroutine calls:                                             CALCP21A.105    
      external calc_rs         ! Compute pseudo radius.                    CALCP21A.106    
      external pv_pint         ! Interpolate variables to theta level.     CALCP21A.107    
      external vortic3         ! Compute term 1.                           CALCP21A.108    
      external vortic4         ! Compute term 2.                           CALCP21A.109    
      external vortic5         ! Compute term 3.                           CALCP21A.110    
      external uv_to_p         ! Interpolate u-grid field to p-grid fld.   CALCP21A.111    
                                                                           CALCP21A.112    
C*----------------------------------------------------------------------   CALCP21A.113    
C*L Call comdecks to get required variables:                               CALCP21A.114    
*CALL C_G                                                                  CALCP21A.115    
*CALL C_OMEGA                                                              CALCP21A.116    
                                                                           CALCP21A.117    
C*----------------------------------------------------------------------   CALCP21A.118    
C*L Define local variables.                                                CALCP21A.119    
      integer i,j,k          ! Loop variables.                             CALCP21A.120    
*IF DEF,MPP                                                                GSM3F403.715    
*IF -DEF,GLOBAL                                                            GSM3F403.716    
      integer START,END                                                    GSM3F403.717    
*ENDIF                                                                     GSM3F403.718    
      INTEGER info !GCG return code                                        GSM3F403.719    
*ENDIF                                                                     GSM3F403.720    
      real mn                ! Mean value used in computing pole values.   CALCP21A.121    
                                                                           CALCP21A.122    
                                                                           CALCP21A.123    
C ----------------------------------------------------------------------   CALCP21A.124    
CL Section 1 Compute rs, the pseudo radius.                                CALCP21A.125    
CL ~~~~~~~~~                                                               CALCP21A.126    
C ----------------------------------------------------------------------   CALCP21A.127    
                                                                           CALCP21A.128    
CL Section 1.1 Call CALC_RS to get rs for level 1.                         CALCP21A.129    
CL ~~~~~~~~~~~ Rs is returned in rs(1,k).                                  CALCP21A.130    
CL             Ts is returned in work1. Rs at level k-1 is input as        CALCP21A.131    
CL             rs(1,2) as at k-1=0 the input is not used by calc_rs.       CALCP21A.132    
                                                                           CALCP21A.133    
      k=1                                                                  CALCP21A.134    
      call calc_rs                                                         CALCP21A.135    
     1            (pstar,ak,bk,work1,rs(1,2),rs(1,k),p_field,k,p_levels,   GSS1F304.154    
     2             llints)                                                 GSS1F304.155    
                                                                           CALCP21A.137    
CL Section 1.2 Call CALC_RS to get rs for level k.                         CALCP21A.138    
CL ~~~~~~~~~~~ Rs is returned in rs(1,k).                                  CALCP21A.139    
CL             Ts is returned in work1. Rs at level k-1 is input as        CALCP21A.140    
CL             rs(1,k-1).                                                  CALCP21A.141    
CL             Loop from 2 to p_levels.                                    CALCP21A.142    
                                                                           CALCP21A.143    
      do 100 k=2,p_levels                                                  CALCP21A.144    
                                                                           CALCP21A.145    
        call calc_rs                                                       CALCP21A.146    
     1          (pstar,ak,bk,work1,rs(1,k-1),rs(1,k),p_field,k,p_levels,   GSS1F304.156    
     2           llints)                                                   GSS1F304.157    
 100  continue                                                             CALCP21A.149    
                                                                           CALCP21A.150    
                                                                           CALCP21A.151    
C ----------------------------------------------------------------------   CALCP21A.152    
CL Section 2 Interpolate p, u, v and rs onto desired theta level,          CALCP21A.153    
CL ~~~~~~~~~ and compute Dtheta/Dp, D(rs.u)/D(p) and D(rs.v)/D(p)          CALCP21A.154    
CL           assuming cubic variation where possible.                      CALCP21A.155    
CL           Interpolate D(rs.u)/D(p) and D(rs.v)/D(p) to the p-grid.      CALCP21A.156    
C ----------------------------------------------------------------------   CALCP21A.157    
                                                                           CALCP21A.158    
      call pv_pint                                                         CALCP21A.159    
     1            (pstar,theta,rs,u,v,p_field,u_field,                     CALCP21A.160    
     2             p_levels,row_length,                                    CALCP21A.161    
*CALL ARGFLDPT                                                             GSM3F403.721    
     3             rmdi,ak,bk,des_press,                                   CALCP21A.162    
     & eta_level,e_levels,dth_dph,n_levels,                                TD141293.180    
     4             theta_on_press,rs_on_press,                             CALCP21A.163    
     5             u_on_press,v_on_press,                                  CALCP21A.164    
     6             drsu_dp,drsv_dp,dtheta_dp)                              CALCP21A.165    
                                                                           CALCP21A.166    
CL Section 2.1 Interpolate D(rs.u)/D(p) and D(rs.v)/D(p) to the p-grid.    CALCP21A.167    
CL ~~~~~~~~~~~                                                             CALCP21A.168    
                                                                           CALCP21A.169    
      do i=1,row_length                                                    GSM3F403.722    
        work1(i)=drsu_dp(i)                                                GSM3F403.723    
      enddo                                                                GSM3F403.724    
      call uv_to_p                                                         CALCP21A.170    
     1            (drsu_dp,work1(row_length+1),u_field,p_field,            GSM3F403.725    
*IF DEF,MPP                                                                GSM3F403.726    
     2             row_length,U_LAST_ROW+1)                                GSM3F403.727    
*ELSE                                                                      GSM3F403.728    
     2             row_length,u_field/row_length)                          CALCP21A.172    
*ENDIF                                                                     GSM3F403.729    
      do i=1,last_p_fld_pt                                                 GSM3F403.730    
        drsu_dp(i)=work1(i)                                                GSM3F403.731    
      enddo                                                                GSM3F403.732    
                                                                           CALCP21A.173    
      do i=1,row_length                                                    GSM3F403.733    
        work1(i)=drsv_dp(i)                                                GSM3F403.734    
      enddo                                                                GSM3F403.735    
      call uv_to_p                                                         CALCP21A.174    
     1            (drsv_dp,work1(row_length+1),u_field,p_field,            GSM3F403.736    
*IF DEF,MPP                                                                GSM3F403.737    
     2             row_length,U_LAST_ROW+1)                                GSM3F403.738    
*ELSE                                                                      GSM3F403.739    
     2             row_length,u_field/row_length)                          CALCP21A.176    
*ENDIF                                                                     GSM3F403.740    
      do i=1,last_p_fld_pt                                                 GSM3F403.741    
        drsv_dp(i)=work1(i)                                                GSM3F403.742    
      enddo                                                                GSM3F403.743    
                                                                           CALCP21A.178    
C ----------------------------------------------------------------------   CALCP21A.179    
CL Section 3 Compute the various terms in the potential vorticity          CALCP21A.180    
CL ~~~~~~~~~ equation for desired theta surface.                           CALCP21A.181    
C ----------------------------------------------------------------------   CALCP21A.182    
                                                                           CALCP21A.183    
CL Section 3.1 Compute 1/(rs*rs*cos(phi))*D(rs.v)/D(p)*D(theta)/D(lamda)   CALCP21A.184    
CL ~~~~~~~~~~~       - 1/(rs*rs)*D(rs.u)/D(p)*D(theta)/D(phi).             CALCP21A.185    
      call vortic3                                                         CALCP21A.186    
     1            (rs_on_press,theta_on_press,drsu_dp,drsv_dp,             CALCP21A.187    
     2             sec_p_latitude,                                         CALCP21A.188    
     3             vorticity3,dtheta_dlatitude,                            CALCP21A.189    
     4             p_field,row_length,                                     CALCP21A.190    
*CALL ARGFLDPT                                                             GSM3F403.744    
     5             latitude_step_inverse,longitude_step_inverse)           CALCP21A.191    
                                                                           CALCP21A.192    
                                                                           CALCP21A.193    
CL Section 3.2 Compute 1/rs*2*omega*cos(phi) * D(theta)/D(phi).            CALCP21A.194    
CL ~~~~~~~~~~~                                                             CALCP21A.195    
      call vortic4                                                         CALCP21A.196    
     1            (rs_on_press,theta_on_press,dtheta_dlatitude,            CALCP21A.197    
     2             sec_p_latitude,                                         CALCP21A.198    
     3             vorticity4,                                             CALCP21A.199    
     4             p_field,row_length,                                     CALCP21A.200    
*CALL ARGFLDPT                                                             GSM3F403.745    
     5             latitude_step_inverse,longitude_step_inverse)           CALCP21A.201    
                                                                           CALCP21A.202    
CL Section 3.3 Compute 1/(rs*cos(phi))*                                    CALCP21A.203    
CL ~~~~~~~~~~~         (D(v)/D(lambda)-D(u.cos(phi))/D(phi)).              CALCP21A.204    
      call vortic5                                                         CALCP21A.205    
     1            (u_on_press,v_on_press,rs_on_press,                      CALCP21A.206    
     2             cos_u_latitude,sec_p_latitude,                          CALCP21A.207    
     3             vorticity5,                                             CALCP21A.208    
     4             p_field,u_field,row_length,                             CALCP21A.209    
*CALL ARGFLDPT                                                             GSM3F403.746    
     5             latitude_step_inverse,longitude_step_inverse)           CALCP21A.210    
                                                                           CALCP21A.211    
                                                                           CALCP21A.212    
C ----------------------------------------------------------------------   CALCP21A.213    
CL Section 4 Compute the potential vorticity using the vorticity terms,    CALCP21A.214    
CL ~~~~~~~~~ D(theta)/D(p) and the value of the coriolis term.             CALCP21A.215    
CL           Care needs to be taken for rmdi.                              CALCP21A.216    
C ----------------------------------------------------------------------   CALCP21A.217    
                                                                           CALCP21A.218    
CL Section 4.1 Firstly, find which points we can actually calculate        CALCP21A.219    
CL ~~~~~~~~~~~ potential vorticity on. To do this we test for missing      CALCP21A.220    
CL             data at all points in a 3x3 stencil around each point.      CALCP21A.221    
CL             Special care is taken at the boundaries and poles.          CALCP21A.222    
                                                                           CALCP21A.223    
C Set logical mask default to .TRUE.                                       CALCP21A.224    
      do i = 1,p_field                                                     CALCP21A.225    
        mask(i) = .true.                                                   CALCP21A.226    
      end do                                                               CALCP21A.227    
                                                                           CALCP21A.228    
C Check points on non-polar rows:                                          CALCP21A.229    
*IF DEF,MPP                                                                GSM3F403.747    
      do j=FIRST_ROW,P_LAST_ROW                                            GSM3F403.748    
*ELSE                                                                      GSM3F403.749    
      do j=2,p_field/row_length-1                                          CALCP21A.230    
*ENDIF                                                                     GSM3F403.750    
C First point on each row:                                                 CALCP21A.231    
        i=(j-1)*row_length+1                                               CALCP21A.232    
*IF DEF,GLOBAL                                                             CALCP21A.233    
*IF -DEF,MPP                                                               GSM3F403.751    
        if (   eta_level(i-1           ).eq.rmdi.                          TD141293.181    
     &      or.eta_level(i-row_length  ).eq.rmdi.                          TD141293.182    
     &      or.eta_level(i-row_length+1).eq.rmdi.                          TD141293.183    
     &      or.eta_level(i+row_length-1).eq.rmdi.                          TD141293.184    
     &      or.eta_level(i             ).eq.rmdi.                          TD141293.185    
     &      or.eta_level(i+1           ).eq.rmdi.                          TD141293.186    
     &      or.eta_level(i+2*row_length-1).eq.rmdi.                        TD141293.187    
     &      or.eta_level(i+row_length  ).eq.rmdi.                          TD141293.188    
     &      or.eta_level(i+row_length+1).eq.rmdi)   mask(i)=.false.        TD141293.189    
*ENDIF                                                                     GSM3F403.752    
*ELSE                                                                      GSM3F403.753    
*IF DEF,MPP                                                                GSM3F403.754    
        IF (at_left_of_LPG) THEN                                           GSM3F403.755    
          mask(i+EW_Halo)=.false.                                          GSM3F403.756    
        ENDIF                                                              GSM3F403.757    
*ELSE                                                                      CALCP21A.243    
        mask(i)=.false.                                                    CALCP21A.244    
*ENDIF                                                                     CALCP21A.245    
*ENDIF                                                                     GSM3F403.758    
C Inner points on each row:                                                CALCP21A.246    
*IF DEF,MPP                                                                GSM3F403.759    
*IF -DEF,GLOBAL                                                            GSM3F403.760    
C Borders of area done separately                                          GSM3F403.761    
        START=EW_Halo+1                                                    GSM3F403.762    
        END=EW_Halo                                                        GSM3F403.763    
        IF (at_left_of_LPG) START=START+1                                  GSM3F403.764    
        IF (at_right_of_LPG) END=END+1                                     GSM3F403.765    
        do i=(j-1)*row_length+START,(j-1)*row_length+row_length-END        GSM3F403.766    
*ELSE                                                                      GSM3F403.767    
        do i=(j-1)*row_length+1+EW_Halo,                                   GSM3F403.768    
     &      (j-1)*row_length+row_length-EW_Halo                            GSM3F403.769    
*ENDIF                                                                     GSM3F403.770    
*ELSE                                                                      GSM3F403.771    
        do i=(j-1)*row_length+2,(j-1)*row_length+row_length-1              CALCP21A.247    
*ENDIF                                                                     GSM3F403.772    
        if (   eta_level(i-row_length-1).eq.rmdi.                          TD141293.190    
     &      or.eta_level(i-row_length  ).eq.rmdi.                          TD141293.191    
     &      or.eta_level(i-row_length+1).eq.rmdi.                          TD141293.192    
     &      or.eta_level(i-1           ).eq.rmdi.                          TD141293.193    
     &      or.eta_level(i             ).eq.rmdi.                          TD141293.194    
     &      or.eta_level(i+1           ).eq.rmdi.                          TD141293.195    
     &      or.eta_level(i+row_length-1).eq.rmdi.                          TD141293.196    
     &      or.eta_level(i+row_length  ).eq.rmdi.                          TD141293.197    
     &      or.eta_level(i+row_length+1).eq.rmdi)   mask(i)=.false.        TD141293.198    
        end do                                                             CALCP21A.257    
C Last point on each row:                                                  CALCP21A.258    
        i=j*row_length                                                     CALCP21A.259    
*IF DEF,GLOBAL                                                             CALCP21A.260    
*IF -DEF,MPP                                                               GSM3F403.773    
        if (   eta_level(i-row_length-1).eq.rmdi.                          TD141293.199    
     &      or.eta_level(i-row_length  ).eq.rmdi.                          TD141293.200    
     &      or.eta_level(i-2*row_length+1).eq.rmdi.                        TD141293.201    
     &      or.eta_level(i-1           ).eq.rmdi.                          TD141293.202    
     &      or.eta_level(i             ).eq.rmdi.                          TD141293.203    
     &      or.eta_level(i-row_length+1).eq.rmdi.                          TD141293.204    
     &      or.eta_level(i+row_length-1).eq.rmdi.                          TD141293.205    
     &      or.eta_level(i+row_length  ).eq.rmdi.                          TD141293.206    
     &      or.eta_level(i+1           ).eq.rmdi)   mask(i)=.false.        TD141293.207    
*ENDIF                                                                     GSM3F403.774    
*ELSE                                                                      GSM3F403.775    
*IF DEF,MPP                                                                GSM3F403.776    
        IF (at_right_of_LPG) THEN                                          GSM3F403.777    
          mask(i-EW_Halo)=.false.                                          GSM3F403.778    
        ENDIF                                                              GSM3F403.779    
*ELSE                                                                      CALCP21A.270    
        mask(i)=.false.                                                    CALCP21A.271    
*ENDIF                                                                     CALCP21A.272    
*ENDIF                                                                     GSM3F403.780    
      end do ! do j = ...                                                  GSM3F403.781    
                                                                           CALCP21A.274    
*IF DEF,GLOBAL                                                             CALCP21A.275    
C Check North pole:                                                        CALCP21A.276    
*IF DEF,MPP                                                                GSM3F403.782    
      IF (at_top_of_LPG) THEN                                              GSM3F403.783    
*ENDIF                                                                     GSM3F403.784    
      if (eta_level(TOP_ROW_START+LAST_ROW_PT-1) .eq. rmdi) then           GSM3F403.785    
        j=1                                                                GSM3F403.786    
      ELSE                                                                 GSM3F403.787    
        do i=TOP_ROW_START+ROW_LENGTH,                                     GSM3F403.788    
     &       TOP_ROW_START+ROW_LENGTH+LAST_ROW_PT-1                        GSM3F403.789    
          if(eta_level(i).eq.rmdi) j=1                                     GSM3F403.790    
        enddo                                                              GSM3F403.791    
      endif                                                                GSM3F403.792    
*IF DEF,MPP                                                                GSM3F403.793    
      CALL GCG_IMAX(1,GC_ROW_GROUP,info,j)                                 GSM3F403.794    
*ENDIF                                                                     GSM3F403.795    
      if (j.eq.1) then                                                     GSM3F403.796    
        do i=TOP_ROW_START,TOP_ROW_START+LAST_ROW_PT-1                     GSM3F403.797    
          mask(i)=.false.                                                  GSM3F403.798    
        enddo                                                              GSM3F403.799    
      endif                                                                GSM3F403.800    
*IF DEF,MPP                                                                GSM3F403.801    
      endif ! if at_top_of_LPG                                             GSM3F403.802    
*ENDIF                                                                     GSM3F403.803    
                                                                           CALCP21A.285    
C Check South pole:                                                        CALCP21A.286    
      j=0                                                                  CALCP21A.287    
*IF DEF,MPP                                                                GSM3F403.804    
      IF (at_base_of_LPG) THEN                                             GSM3F403.805    
*ENDIF                                                                     GSM3F403.806    
      if (eta_level(P_BOT_ROW_START) .eq. rmdi) then                       GSM3F403.807    
        j=1                                                                GSM3F403.808    
      else                                                                 GSM3F403.809    
        do i=P_BOT_ROW_START-ROW_LENGTH,                                   GSM3F403.810    
     &       P_BOT_ROW_START-ROW_LENGTH+LAST_ROW_PT-1                      GSM3F403.811    
          if (eta_level(i) .eq. rmdi) j=1                                  GSM3F403.812    
        enddo                                                              GSM3F403.813    
      endif                                                                GSM3F403.814    
*IF DEF,MPP                                                                GSM3F403.815    
      CALL GCG_IMAX(1,GC_ROW_GROUP,info,j)                                 GSM3F403.816    
*ENDIF                                                                     GSM3F403.817    
                                                                           GSM3F403.818    
      if (j.eq.1) then                                                     GSM3F403.819    
        do i=P_BOT_ROW_START,P_BOT_ROW_START+LAST_ROW_PT-1                 GSM3F403.820    
          mask(i)=.false.                                                  GSM3F403.821    
        enddo                                                              GSM3F403.822    
      endif                                                                GSM3F403.823    
*IF DEF,MPP                                                                GSM3F403.824    
      endif ! if at_base_of_LPG                                            GSM3F403.825    
*ENDIF                                                                     GSM3F403.826    
*ELSE                                                                      CALCP21A.295    
*IF DEF,MPP                                                                GSM3F403.827    
      IF (at_top_of_LPG) THEN                                              GSM3F403.828    
*ENDIF                                                                     GSM3F403.829    
        do i=TOP_ROW_START,TOP_ROW_START+LAST_ROW_PT-1                     GSM3F403.830    
          mask(i)=.false.                                                  GSM3F403.831    
        enddo                                                              GSM3F403.832    
*IF DEF,MPP                                                                GSM3F403.833    
      ENDIF ! if at_top_of_LPG                                             GSM3F403.834    
                                                                           GSM3F403.835    
      IF (at_base_of_LPG) THEN                                             GSM3F403.836    
*ENDIF                                                                     GSM3F403.837    
        do i=P_BOT_ROW_START,P_BOT_ROW_START+LAST_ROW_PT-1                 GSM3F403.838    
          mask(i)=.false.                                                  GSM3F403.839    
        enddo                                                              GSM3F403.840    
*IF DEF,MPP                                                                GSM3F403.841    
      ENDIF ! if at_base_of_LPG                                            GSM3F403.842    
*ENDIF                                                                     GSM3F403.843    
*ENDIF                                                                     CALCP21A.303    
                                                                           CALCP21A.304    
                                                                           CALCP21A.305    
CL Section 4.2 Interpolate F3 (the Coriolis term) to the p-grid.           CALCP21A.306    
CL ~~~~~~~~~~~                                                             CALCP21A.307    
                                                                           CALCP21A.308    
      call uv_to_p                                                         CALCP21A.309    
     1            (f3,f3_p(row_length+1),u_field,p_field,                  CALCP21A.310    
*IF DEF,MPP                                                                GSM3F403.844    
     2             row_length,U_LAST_ROW+1)                                GSM3F403.845    
*ELSE                                                                      GSM3F403.846    
     2             row_length,u_field/row_length)                          CALCP21A.311    
*ENDIF                                                                     GSM3F403.847    
                                                                           CALCP21A.312    
                                                                           CALCP21A.313    
CL Section 4.3 Calculate the potential vorticity.                          CALCP21A.314    
CL ~~~~~~~~~~~                                                             CALCP21A.315    
                                                                           CALCP21A.316    
      do i = START_POINT_NO_HALO,END_P_POINT_NO_HALO                       GSM3F403.848    
        if (mask(i)) then                                                  CALCP21A.318    
          pvort_p(i) = (vorticity3(i) + vorticity4(i) -                    CALCP21A.319    
     &                  dtheta_dp(i)*(f3_p(i) + vorticity5(i))) * g        CALCP21A.320    
        else                                                               CALCP21A.321    
          pvort_p(i) = rmdi                                                CALCP21A.322    
        end if                                                             CALCP21A.323    
      enddo                                                                CALCP21A.324    
                                                                           CALCP21A.325    
CL Section 4.3.1 Compute the potential vorticity at the Northern           CALCP21A.326    
CL ~~~~~~~~~~~~~ and Southern boundaries.                                  CALCP21A.327    
CL               Note that under UPDATE identifier GLOBAL, as pv           CALCP21A.328    
CL               is a scaler this is defined to be the mean of             CALCP21A.329    
CL               the near pole row.                                        CALCP21A.330    
CL               If any rmdi is present in the row next to the pole,       CALCP21A.331    
CL               then the pole has a value rmdi.                           CALCP21A.332    
CL               ELSE, the value returned is missing data.                 CALCP21A.333    
                                                                           CALCP21A.334    
C Work out the Northern boundary.                                          CALCP21A.335    
*IF DEF,MPP                                                                GSM3F403.849    
      if (at_top_of_LPG) then                                              GSM3F403.850    
*ENDIF                                                                     GSM3F403.851    
        if (mask(TOP_ROW_START+LAST_ROW_PT-1)) then                        GSM3F403.852    
                                                                           CALCP21A.352    
          mn=0                                                             GSM3F403.853    
*IF -DEF,MPP                                                               GSM3F403.854    
          do i=TOP_ROW_START+row_length,                                   GSM3F403.855    
     &         TOP_ROW_START+2*row_length-1                                GSM3F403.856    
            mn=mn+pvort_p(i)                                               GSM3F403.857    
          enddo                                                            GSM3F403.858    
*ELSE                                                                      GSM3F403.859    
*IF DEF,REPROD                                                             GSM3F403.860    
          call GCG_RVECSUMR(                                               GSM3F403.861    
*ELSE                                                                      GSM3F403.862    
          call GCG_RVECSUMF(                                               GSM3F403.863    
*ENDIF                                                                     GSM3F403.864    
     &      row_length-2*EW_Halo,row_length-2*EW_Halo,1,1,                 GSM3F403.865    
     &      pvort_p(TOP_ROW_START+FIRST_ROW_PT-1+row_length),              GSM3F403.866    
     &      GC_ROW_GROUP,info,mn)                                          GSM3F403.867    
*ENDIF                                                                     GSM3F403.868    
          mn=mn/GLOBAL_ROW_LENGTH                                          GSM3F403.869    
                                                                           CALCP21A.370    
          do i=TOP_ROW_START,TOP_ROW_START+row_length-1                    GSM3F403.870    
            pvort_p(i)=mn                                                  GSM3F403.871    
          enddo                                                            GSM3F403.872    
        else                                                               GSM3F403.873    
! Missing data                                                             GSM3F403.874    
          do i=TOP_ROW_START,TOP_ROW_START+row_length-1                    GSM3F403.875    
            pvort_p(i)=rmdi                                                GSM3F403.876    
          enddo                                                            GSM3F403.877    
        endif                                                              GSM3F403.878    
*IF DEF,MPP                                                                GSM3F403.879    
      endif ! if at_top_of_LPG                                             GSM3F403.880    
*ENDIF                                                                     GSM3F403.881    
                                                                           GSM3F403.882    
! Southern boundary                                                        GSM3F403.883    
*IF DEF,MPP                                                                GSM3F403.884    
                                                                           GSM3F403.885    
      if (at_base_of_LPG) then                                             GSM3F403.886    
*ENDIF                                                                     GSM3F403.887    
        if (mask(P_BOT_ROW_START+FIRST_ROW_PT-1)) then                     GSM3F403.888    
                                                                           GSM3F403.889    
          mn=0                                                             GSM3F403.890    
*IF -DEF,MPP                                                               GSM3F403.891    
          do i=P_BOT_ROW_START-ROW_LENGTH,                                 GSM3F403.892    
     &         P_BOT_ROW_START-1                                           GSM3F403.893    
            mn=mn+pvort_p(i)                                               GSM3F403.894    
          enddo                                                            GSM3F403.895    
*ELSE                                                                      GSM3F403.896    
*IF DEF,REPROD                                                             GSM3F403.897    
          call GCG_RVECSUMR(                                               GSM3F403.898    
*ELSE                                                                      GSM3F403.899    
          call GCG_RVECSUMF(                                               GSM3F403.900    
*ENDIF                                                                     GSM3F403.901    
     &      row_length-2*EW_Halo,row_length-2*EW_Halo,1,1,                 GSM3F403.902    
     &      pvort_p(P_BOT_ROW_START+FIRST_ROW_PT-1-row_length),            GSM3F403.903    
     &      GC_ROW_GROUP,info,mn)                                          GSM3F403.904    
*ENDIF                                                                     GSM3F403.905    
          mn=mn/GLOBAL_ROW_LENGTH                                          GSM3F403.906    
                                                                           GSM3F403.907    
          do i=P_BOT_ROW_START,P_BOT_ROW_START+row_length-1                GSM3F403.908    
            pvort_p(i)=mn                                                  GSM3F403.909    
          enddo                                                            GSM3F403.910    
        else                                                               GSM3F403.911    
! Missing data                                                             GSM3F403.912    
          do i=P_BOT_ROW_START,P_BOT_ROW_START+row_length-1                GSM3F403.913    
            pvort_p(i)=rmdi                                                GSM3F403.914    
          enddo                                                            GSM3F403.915    
        endif                                                              GSM3F403.916    
*IF DEF,MPP                                                                GSM3F403.917    
      endif ! if at_base_of_LPG                                            GSM3F403.918    
*ENDIF                                                                     GSM3F403.919    
                                                                           GSM3F403.920    
*IF DEF,MPP                                                                GSM3F403.921    
! Fill in the N+S halo areas with rmdi                                     GSM3F403.922    
      do i=1,FIRST_FLD_PT-1                                                GSM3F403.923    
        pvort_p(i)=rmdi                                                    GSM3F403.924    
      enddo                                                                GSM3F403.925    
      do i=LAST_P_FLD_PT+1,p_field                                         GSM3F403.926    
        pvort_p(i)=rmdi                                                    GSM3F403.927    
      enddo                                                                GSM3F403.928    
*ENDIF                                                                     GSM3F403.929    
                                                                           CALCP21A.371    
      return                                                               CALCP21A.372    
      end                                                                  CALCP21A.373    
                                                                           CALCP21A.374    
*ENDIF                                                                     CALCP21A.375