*IF DEF,A12_1C,AND,DEF,MPP                                                 THADV1C.2      
C *****************************COPYRIGHT******************************     THADV1C.3      
C (c) CROWN COPYRIGHT 1996, METEOROLOGICAL OFFICE, All Rights Reserved.    THADV1C.4      
C                                                                          THADV1C.5      
C Use, duplication or disclosure of this code is subject to the            THADV1C.6      
C restrictions as set forth in the contract.                               THADV1C.7      
C                                                                          THADV1C.8      
C                Meteorological Office                                     THADV1C.9      
C                London Road                                               THADV1C.10     
C                BRACKNELL                                                 THADV1C.11     
C                Berkshire UK                                              THADV1C.12     
C                RG12 2SZ                                                  THADV1C.13     
C                                                                          THADV1C.14     
C If no contract has been raised with this copy of the code, the use,      THADV1C.15     
C duplication or disclosure of it is strictly prohibited.  Permission      THADV1C.16     
C to do so must first be obtained in writing from the Head of Numerical    THADV1C.17     
C Modelling at the above address.                                          THADV1C.18     
C ******************************COPYRIGHT******************************    THADV1C.19     
CLL   SUBROUTINE TH_ADV -------------------------------------------        THADV1C.20     
CLL                                                                        THADV1C.21     
CLL   PURPOSE:  CALCULATES MASS-WEIGHTED INCREMENTS TO THETAL              THADV1C.22     
CLL DUE TO ADVECTION  BY USING EQUATION (35) TO CALCULATE PROVISIONAL      THADV1C.23     
CLL VALUES OF THETAL AT THE NEW TIME-LEVEL, AND THEN RECALCULATING THE     THADV1C.24     
CLL ADVECTION TERMS ON THE RIGHT-HAND SIDE OF (35) USING THESE             THADV1C.25     
CLL PROVISIONAL VALUES. THE FINAL INCREMENTS ARE CALCULATED AS IN          THADV1C.26     
CLL EQUATION (40). THOSE REQUIRING FILTERING ARE FILTERED AND ALL THE      THADV1C.27     
CLL INCREMENTS ARE ADDED ONTO THE FIELDS USING (40).  IF RUNNING A         THADV1C.28     
CLL GLOBAL MODEL POLAR IS CALLED TO UPDATE POLAR VALUES.                   THADV1C.29     
CLL                                                                        THADV1C.30     
CLL   NOT SUITABLE FOR SINGLE COLUMN USE.                                  THADV1C.31     
CLL   VERSION FOR CRAY Y-MP                                                THADV1C.32     
CLL                                                                        THADV1C.33     
CLL MM, DR      <- PROGRAMMER OF SOME OR ALL OF PREVIOUS CODE OR CHANGES   THADV1C.34     
CLL MPP CODE ADDED BY P.BURTON                                             THADV1C.35     
CLL                                                                        THADV1C.36     
CLL  MODEL            MODIFICATION HISTORY FROM MODEL VERSION 4.1:         THADV1C.37     
CLL VERSION  DATE                                                          THADV1C.38     
CLL                                                                        THADV1C.39     
CLL 4.1      07/12/95 New version of routine specifically for MPP          THADV1C.40     
CLL                   P.Burton                                             THADV1C.41     
!LL   4.2    16/08/96  Add TYPFLDPT arguments to FILTER subroutine         APB0F402.42     
!LL                    and make the FILTER_WAVE_NUMBER arrays              APB0F402.43     
!LL                    globally sized                       P.Burton       APB0F402.44     
!LL   4.2    10/01/97  Initialise unprocessed points in THETAL_PROV.       ADR2F402.1      
!LL                    D. Robinson.                                        ADR2F402.2      
!LL 4.3      24/04/97 Fixes to 4th order calculations   P.Burton           GPB5F403.17     
C     vn4.3    Mar. 97   T3E migration : optimisation changes              GSS1F403.706    
C                                       D.Salmond                          GSS1F403.707    
CLL                                                                        THADV1C.42     
CLL   PROGRAMMING STANDARD: UNIFIED MODEL DOCUMENTATION PAPER NO. 4,       THADV1C.43     
CLL                         STANDARD B.                                    THADV1C.44     
CLL                                                                        THADV1C.45     
CLL   SYSTEM COMPONENTS COVERED: P121                                      THADV1C.46     
CLL                                                                        THADV1C.47     
CLL   SYSTEM TASK: P1                                                      THADV1C.48     
CLL                                                                        THADV1C.49     
CLL   DOCUMENTATION:       THE EQUATIONS USED ARE (35) AND (40)            THADV1C.50     
CLL                        IN UNIFIED MODEL DOCUMENTATION PAPER NO. 10     THADV1C.51     
CLL                        M.J.P. CULLEN,T.DAVIES AND M.H.MAWSON           THADV1C.52     
CLL                                                                        THADV1C.53     
CLLEND-------------------------------------------------------------        THADV1C.54     
                                                                           THADV1C.55     
C*L   ARGUMENTS:---------------------------------------------------        THADV1C.56     

      SUBROUTINE TH_ADV                                                     1,59THADV1C.57     
     1              (THETAL,PSTAR_OLD,PSTAR,U_MEAN,V_MEAN,                 THADV1C.58     
     2              SEC_P_LATITUDE,ETADOT_MEAN,RS,DELTA_AK,DELTA_BK,       THADV1C.59     
     3              LATITUDE_STEP_INVERSE,ADVECTION_TIMESTEP,NU_BASIC,     THADV1C.60     
     4              LONGITUDE_STEP_INVERSE,NORTHERN_FILTERED_P_ROW,        THADV1C.61     
     5              SOUTHERN_FILTERED_P_ROW,P_LEVELS,                      THADV1C.62     
     6              U_FIELD,P_FIELD,ROW_LENGTH,                            THADV1C.63     
*CALL ARGFLDPT                                                             THADV1C.64     
     7              TRIGS,IFAX,FILTER_WAVE_NUMBER_P_ROWS,SEC_U_LATITUDE,   THADV1C.65     
     8              AKH,BKH,QCL,QCF,P_EXNER,OMEGA,                         THADV1C.66     
     9              Q_LEVELS,AK,BK,L_SECOND,                               THADV1C.67     
     &              extended_address,                                      GSS1F403.708    
     &              LWHITBROM)                                             THADV1C.68     
                                                                           THADV1C.69     
      IMPLICIT NONE                                                        THADV1C.70     
                                                                           THADV1C.71     
! All TYPFLDPT arguments are intent IN                                     THADV1C.72     
*CALL TYPFLDPT                                                             THADV1C.73     
                                                                           THADV1C.74     
      INTEGER                                                              THADV1C.75     
     *  P_FIELD            !IN DIMENSION OF FIELDS ON PRESSSURE GRID.      THADV1C.76     
     *, U_FIELD            !IN DIMENSION OF FIELDS ON VELOCITY GRID        THADV1C.77     
     *, P_LEVELS           !IN NUMBER OF PRESSURE LEVELS.                  THADV1C.78     
     *, Q_LEVELS           !IN NUMBER OF MOIST LEVELS.                     THADV1C.79     
     *, ROW_LENGTH         !IN NUMBER OF POINTS PER ROW                    THADV1C.80     
     *, NORTHERN_FILTERED_P_ROW !IN ROW ON WHICH FILTERING STOPS           THADV1C.81     
     *, SOUTHERN_FILTERED_P_ROW !IN ROW ON WHICH FILTERING STARTS AGAIN.   THADV1C.82     
     &, FILTER_WAVE_NUMBER_P_ROWS(GLOBAL_P_FIELD/GLOBAL_ROW_LENGTH)        APB0F402.45     
     &       ! LAST WAVE NUMBER NOT TO BE CHOPPED                          APB0F402.46     
     *, IFAX(10)           !IN HOLDS FACTORS OF ROW_LENGTH USED BY         THADV1C.85     
     *                     ! FILTERING.                                    THADV1C.86     
                                                                           THADV1C.87     
C LOGICAL VARIABLE                                                         THADV1C.88     
      LOGICAL                                                              THADV1C.89     
     *  L_SECOND     ! SET TO TRUE IF NU_BASIC IS ZERO.                    THADV1C.90     
     & ,LWHITBROM                                                          THADV1C.91     
      INTEGER extended_address(P_FIELD)                                    GSS1F403.709    
                                                                           THADV1C.92     
      REAL                                                                 THADV1C.93     
     * THETAL(P_FIELD,P_LEVELS)  !INOUT THETAL FIELD                       THADV1C.94     
     *                           ! MASS-WEIGHTED ON OUTPUT.                THADV1C.95     
                                                                           THADV1C.96     
      REAL                                                                 THADV1C.97     
     * U_MEAN(U_FIELD,P_LEVELS) !IN AVERAGED MASS-WEIGHTED U VELOCITY      THADV1C.98     
     *                          !   FROM ADJUSTMENT STEP                   THADV1C.99     
     *,V_MEAN(U_FIELD,P_LEVELS) !IN AVERAGED MASS-WEIGHTED V VELOCITY      THADV1C.100    
     *                          !   * COS(LAT) FROM ADJUSTMENT STEP        THADV1C.101    
     *,ETADOT_MEAN(P_FIELD,P_LEVELS)  !IN AVERAGED MASS-WEIGHTED           THADV1C.102    
     *                          !VERTICAL VELOCITY FROM ADJUSTMENT STEP    THADV1C.103    
     *,PSTAR(P_FIELD)           !IN PSTAR FIELD AT NEW TIME-LEVEL          THADV1C.104    
     *,PSTAR_OLD(P_FIELD)       !IN PSTAR AT PREVIOUS TIME-LEVEL           THADV1C.105    
     *,RS(P_FIELD,P_LEVELS)     !IN RS FIELD                               THADV1C.106    
     *,QCL(P_FIELD,Q_LEVELS)    !IN. PRIMARY ARRAY FOR QCL                 THADV1C.107    
     *,QCF(P_FIELD,Q_LEVELS)    !IN. PRIMARY ARRAY FOR QCF                 THADV1C.108    
     *,OMEGA(U_FIELD,P_LEVELS)  !IN. TRUE VERTICAL VELOCITY DP/DT          THADV1C.109    
     *,P_EXNER(P_FIELD,P_LEVELS+1) !IN. PRIMARY ARRAY FOR EXNER FUNCTION   THADV1C.110    
                                                                           THADV1C.111    
      REAL                                                                 THADV1C.112    
     * DELTA_AK(P_LEVELS)      !IN LAYER THICKNESS                         THADV1C.113    
     *,DELTA_BK(P_LEVELS)      !IN LAYER THICKNESS                         THADV1C.114    
     *,AK(P_LEVELS)            !IN HYBRID CO-ORDINATE AT FULL LEVELS       THADV1C.115    
     *,BK(P_LEVELS)            !IN HYBRID CO-ORDINATE AT FULL LEVELS       THADV1C.116    
     *,AKH(P_LEVELS+1)         !IN HYBRID CO-ORDINATE AT HALF LEVELS       THADV1C.117    
     *,BKH(P_LEVELS+1)         !IN HYBRID CO-ORDINATE AT HALF LEVELS       THADV1C.118    
     *,SEC_P_LATITUDE(P_FIELD) !IN 1/COS(LAT) AT P POINTS (2-D ARRAY)      THADV1C.119    
     *,SEC_U_LATITUDE(U_FIELD) !IN 1/COS(LAT) AT U POINTS (2-D ARRAY)      THADV1C.120    
     *,LONGITUDE_STEP_INVERSE  !IN 1/(DELTA LAMDA)                         THADV1C.121    
     *,LATITUDE_STEP_INVERSE   !IN 1/(DELTA PHI)                           THADV1C.122    
     *,ADVECTION_TIMESTEP      !IN                                         THADV1C.123    
     *,NU_BASIC                !IN STANDARD NU TERM FOR MODEL RUN.         THADV1C.124    
     *,TRIGS(ROW_LENGTH)       !IN HOLDS TRIGONOMETRIC FUNCTIONS USED      THADV1C.125    
     *                         ! IN FILTERING.                             THADV1C.126    
                                                                           THADV1C.127    
C*L   DEFINE ARRAYS AND VARIABLES USED IN THIS ROUTINE-----------------    THADV1C.128    
C DEFINE LOCAL ARRAYS: 24 ARE REQUIRED.                                    THADV1C.129    
      REAL                                                                 THADV1C.130    
     &    OMEGA_P(P_FIELD,P_LEVELS)    ! HOLDS OMEGA AT P POINTS.          THADV1C.131    
                                                                           THADV1C.132    
      REAL                                                                 THADV1C.133    
     * THETAL_FIRST_INC(P_FIELD,P_LEVELS) ! HOLDS THETAL INCREMENT         THADV1C.134    
     *                           ! RETURNED BY FIRST CALL TO ADV_P_GD      THADV1C.135    
     *,THETAL_SECOND_INC(P_FIELD)! HOLDS THETAL INCREMENT                  THADV1C.136    
     *                           ! RETURNED BY SECOND CALL TO ADV_P_GD     THADV1C.137    
     *,THETAL_PROV(P_FIELD,P_LEVELS) ! HOLDS PROVISIONAL VALUE OF          THADV1C.138    
     *                           ! THETAL                                  THADV1C.139    
                                                                           THADV1C.140    
      REAL                                                                 THADV1C.141    
     * THETAL_INCREMENT(P_FIELD,P_LEVELS) !HOLDS INCREMENT TO THETAL       THADV1C.142    
     *,ZERO(P_FIELD)             !A FIELD OF ZEROES USED WHERE VERTICAL    THADV1C.143    
     *                           !VELOCITY IS ZERO.                        THADV1C.144    
                                                                           THADV1C.145    
      REAL                                                                 THADV1C.146    
     * NUX(P_FIELD,P_LEVELS) !COURANT NBR DEPENDENT NU AT P PTS USED       THADV1C.147    
     *                    ! IN EAST-WEST ADVECTION.                        THADV1C.148    
     *,NUY(P_FIELD,P_LEVELS) !COURANT NBR DEPENDENT NU AT P PTS USED       THADV1C.149    
     *                    ! IN NORTH-SOUTH ADVECTION.                      THADV1C.150    
                                                                           THADV1C.151    
      REAL NUX_MIN(upd_P_ROWS), ! minimum value of NUX along a row         THADV1C.152    
     &     NUY_MIN(ROW_LENGTH)  ! min of NUY along a column                THADV1C.153    
                                                                           THADV1C.154    
      REAL                                                                 THADV1C.155    
     * BRSP(P_FIELD,P_LEVELS) !MASS TERM AT LEVEL K                        THADV1C.156    
                                                                           THADV1C.157    
! Work space required to allow the use of Fourth Order Advection           THADV1C.158    
! U/V_MEAN_COPY and T_COPY arrays are defined with an extra halo           THADV1C.159    
! this is required for the bigger stencil of the 4th order operator.       THADV1C.160    
      REAL U_MEAN_COPY((ROW_LENGTH+2*extra_EW_Halo)*                       THADV1C.161    
     &                 (tot_U_ROWS+2*extra_NS_Halo),P_LEVELS),             THADV1C.162    
     &  !    Copy of U_MEAN with extra halo space for 4th order            THADV1C.163    
     &      V_MEAN_COPY((ROW_LENGTH+2*extra_EW_Halo)*                      THADV1C.164    
     &                  (tot_U_ROWS+2*extra_NS_Halo),P_LEVELS),            THADV1C.165    
     &  !    Copy of V_MEAN with extra halo space for 4th order            THADV1C.166    
     &      T_COPY((ROW_LENGTH+2*extra_EW_Halo)*                           THADV1C.167    
     &             (tot_P_ROWS+2*extra_NS_Halo),P_LEVELS)                  THADV1C.168    
     &  !    Copy of THETAL with extra halo space for 4th order            THADV1C.169    
                                                                           THADV1C.170    
      INTEGER  extended_P_FIELD,                                           THADV1C.171    
     &         extended_U_FIELD                                            THADV1C.172    
!  These are the sizes of the arrays with the extra halos                  THADV1C.173    
C*---------------------------------------------------------------------    THADV1C.174    
C DEFINE LOCAL VARIABLES                                                   THADV1C.175    
      INTEGER                                                              THADV1C.176    
     *  P_POINTS_UPDATE    ! NUMBER OF P POINTS TO BE UPDATED.             THADV1C.177    
     *                     !  = ROWS*ROWLENGTH                             THADV1C.178    
     *, U_POINTS_UPDATE    ! NUMBER OF U POINTS TO BE UPDATED.             THADV1C.179    
     *                     !  = (ROWS-1)*ROWLENGTH                         THADV1C.180    
     *, P_POINTS_REQUIRED  ! NUMBER OF P POINTS AT WHICH VALUES ARE        THADV1C.181    
     *                     ! NEEDED TO UPDATE AT P_POINTS_UPDATE           THADV1C.182    
     *, U_POINTS_REQUIRED  ! NUMBER OF U POINTS AT WHICH VALUES ARE        THADV1C.183    
     *                     ! NEEDED TO UPDATE AT U_POINTS_UPDATE           THADV1C.184    
     *, START_U_REQUIRED   ! FIRST U POINT OF VALUES REQUIRED TO UPDATE    THADV1C.185    
     *                     ! AT P POINTS UPDATE.                           THADV1C.186    
     *, END_U_REQUIRED     ! LAST U POINT OF REQUIRED VALUES.              THADV1C.187    
                                                                           THADV1C.188    
      INTEGER I_start,I_end  ! loop bounds                                 THADV1C.189    
      INTEGER info  ! return code from comms                               THADV1C.190    
                                                                           THADV1C.191    
C REAL SCALARS                                                             THADV1C.192    
      REAL                                                                 THADV1C.193    
     & SCALAR1,SCALAR2,CONST1,LC_LF,TIMESTEP                               THADV1C.194    
     &,PK, PK1         ! Pressure at half levels k and k1 (k1=k-1)         THADV1C.195    
     &,P_EXNER_FULL    ! Exner pressure at full model level                THADV1C.196    
                                                                           THADV1C.197    
C COUNT VARIABLES FOR DO LOOPS ETC.                                        THADV1C.198    
      INTEGER                                                              THADV1C.199    
     &  I,J,K1,IK,K                                                        THADV1C.200    
     *, FILTER_SPACE ! HORIZONTAL DIMENSION OF SPACE NEEDED IN FILTERING   THADV1C.201    
     *               ! ROUTINE.                                            THADV1C.202    
                                                                           THADV1C.203    
C*L   EXTERNAL SUBROUTINE CALLS:---------------------------------------    THADV1C.204    
      EXTERNAL ADV_P_GD,POLAR,UV_TO_P,FILTER                               THADV1C.205    
*IF DEF,CRAY                                                               THADV1C.206    
      INTEGER ISMIN                                                        THADV1C.207    
      EXTERNAL ISMIN                                                       THADV1C.208    
*ENDIF                                                                     THADV1C.209    
C*---------------------------------------------------------------------    THADV1C.210    
CL    CALL COMDECK TO GET PHYSICAL CONSTANTS USED.                         THADV1C.211    
                                                                           THADV1C.212    
*CALL C_THADV                                                              THADV1C.213    
                                                                           THADV1C.214    
CL    MAXIMUM VECTOR LENGTH ASSUMED IS P_FIELD.                            THADV1C.215    
CL---------------------------------------------------------------------    THADV1C.216    
CL    INTERNAL STRUCTURE INCLUDING SUBROUTINE CALLS:                       THADV1C.217    
CL---------------------------------------------------------------------    THADV1C.218    
CL                                                                         THADV1C.219    
                                                                           THADV1C.220    
*CALL P_EXNERC                                                             THADV1C.221    
                                                                           THADV1C.222    
CL---------------------------------------------------------------------    THADV1C.223    
CL    SECTION 1.     INITIALISATION                                        THADV1C.224    
CL---------------------------------------------------------------------    THADV1C.225    
C INCLUDE LOCAL CONSTANTS FROM GENERAL CONSTANTS BLOCK                     THADV1C.226    
                                                                           THADV1C.227    
      LC_LF = LC + LF                                                      THADV1C.228    
      P_POINTS_UPDATE   = upd_P_ROWS*ROW_LENGTH                            THADV1C.229    
      U_POINTS_UPDATE   = upd_U_ROWS*ROW_LENGTH                            THADV1C.230    
      P_POINTS_REQUIRED = (upd_P_ROWS+2)*ROW_LENGTH                        THADV1C.231    
      U_POINTS_REQUIRED = (upd_U_ROWS+2)*ROW_LENGTH                        THADV1C.232    
      START_U_REQUIRED  = START_POINT_NO_HALO-ROW_LENGTH                   THADV1C.233    
      END_U_REQUIRED    = END_U_POINT_NO_HALO+ROW_LENGTH                   THADV1C.234    
                                                                           THADV1C.235    
                                                                           THADV1C.236    
C *IF -DEF,NOWHBR replaced by LWHITBROM logical                            THADV1C.237    
      IF (LWHITBROM) THEN                                                  THADV1C.238    
CL    CALCULATE BRSP TERM AT LEVEL K                                       THADV1C.239    
                                                                           THADV1C.240    
      K=1                                                                  THADV1C.241    
! Loop over entire field                                                   THADV1C.242    
      DO I=FIRST_VALID_PT,LAST_P_VALID_PT                                  THADV1C.243    
        BRSP(I,K)=(3.*RS(I,K)+RS(I,K+1))*(RS(I,K)-RS(I,K+1))               THADV1C.244    
     *                *BKH(K+1)*.25*(PSTAR(I)-PSTAR_OLD(I))                THADV1C.245    
      ENDDO                                                                THADV1C.246    
      K=P_LEVELS                                                           THADV1C.247    
! Loop over entire field                                                   THADV1C.248    
      DO I=FIRST_VALID_PT,LAST_P_VALID_PT                                  THADV1C.249    
        BRSP(I,K)=-(3.*RS(I,K)+RS(I,K-1))*(RS(I,K)-RS(I,K-1))              THADV1C.250    
     *                *BKH(K)*.25*(PSTAR(I)-PSTAR_OLD(I))                  THADV1C.251    
      ENDDO                                                                THADV1C.252    
                                                                           THADV1C.253    
      DO K=2,P_LEVELS -1                                                   THADV1C.254    
! Loop over entire field                                                   THADV1C.255    
        DO I=FIRST_VALID_PT,LAST_P_VALID_PT                                THADV1C.256    
          BRSP(I,K)=((3.*RS(I,K)+RS(I,K+1))*(RS(I,K)-RS(I,K+1))*BKH(K+1)   THADV1C.257    
     *              *.25*(PSTAR(I)-PSTAR_OLD(I)))                          THADV1C.258    
     *              -((3.*RS(I,K)+RS(I,K-1))*(RS(I,K)-RS(I,K-1))*BKH(K)    THADV1C.259    
     *              *.25*(PSTAR(I)-PSTAR_OLD(I)))                          THADV1C.260    
        ENDDO                                                              THADV1C.261    
                                                                           THADV1C.262    
      ENDDO                                                                THADV1C.263    
      END IF                                                               THADV1C.264    
C *ENDIF                                                                   THADV1C.265    
                                                                           THADV1C.266    
      DO I=FIRST_VALID_PT,LAST_P_VALID_PT                                  THADV1C.267    
        ZERO(I) = 0.                                                       THADV1C.268    
      ENDDO                                                                THADV1C.269    
                                                                           THADV1C.270    
! In order to use the same call to adv_p_gd for both the second and        THADV1C.271    
! fourth order advection, U/V_MEAN are copied into _COPY arrays.           THADV1C.272    
! In the case of second order advection some of the work space is          THADV1C.273    
! wasted as there is more halo than we need.                               THADV1C.274    
                                                                           THADV1C.275    
! Calculate the size of the extended arrays which contain an               THADV1C.276    
! extra halo:                                                              THADV1C.277    
      extended_U_FIELD=(ROW_LENGTH+2*extra_EW_Halo)*                       THADV1C.278    
     &                 (tot_U_ROWS+2*extra_NS_Halo)                        THADV1C.279    
      extended_P_FIELD=(ROW_LENGTH+2*extra_EW_Halo)*                       THADV1C.280    
     &                 (tot_P_ROWS+2*extra_NS_Halo)                        THADV1C.281    
                                                                           THADV1C.282    
      IF (L_SECOND) THEN                                                   THADV1C.283    
                                                                           THADV1C.284    
! Copy U/V_MEAN to U/V_MEAN_COPY with the same sized halos                 THADV1C.285    
        CALL COPY_FIELD(U_MEAN,U_MEAN_COPY,                                THADV1C.286    
     &                  U_FIELD,extended_U_FIELD,                          THADV1C.287    
     &                  ROW_LENGTH,tot_U_ROWS,P_LEVELS,                    THADV1C.288    
     &                  EW_Halo,NS_Halo,                                   THADV1C.289    
     &                  EW_Halo,NS_Halo,                                   THADV1C.290    
     &                  .FALSE.)                                           THADV1C.291    
        CALL COPY_FIELD(V_MEAN,V_MEAN_COPY,                                THADV1C.292    
     &                  U_FIELD,extended_U_FIELD,                          THADV1C.293    
     &                  ROW_LENGTH,tot_U_ROWS,P_LEVELS,                    THADV1C.294    
     &                  EW_Halo,NS_Halo,                                   THADV1C.295    
     &                  EW_Halo,NS_Halo,                                   THADV1C.296    
     &                  .FALSE.)                                           THADV1C.297    
                                                                           THADV1C.298    
      ELSE  ! if its fourth order:                                         THADV1C.299    
                                                                           THADV1C.300    
        CALL COPY_FIELD(U_MEAN,U_MEAN_COPY,                                THADV1C.301    
     &                  U_FIELD,extended_U_FIELD,                          THADV1C.302    
     &                  ROW_LENGTH,tot_U_ROWS,P_LEVELS,                    THADV1C.303    
     &                  EW_Halo,NS_Halo,                                   THADV1C.304    
     &                  halo_4th,halo_4th,                                 THADV1C.305    
     &                  .TRUE.)                                            THADV1C.306    
        CALL COPY_FIELD(V_MEAN,V_MEAN_COPY,                                THADV1C.307    
     &                  U_FIELD,extended_U_FIELD,                          THADV1C.308    
     &                  ROW_LENGTH,tot_U_ROWS,P_LEVELS,                    THADV1C.309    
     &                  EW_Halo,NS_Halo,                                   THADV1C.310    
     &                  halo_4th,halo_4th,                                 THADV1C.311    
     &                  .TRUE.)                                            THADV1C.312    
        CALL COPY_FIELD(THETAL,T_COPY,                                     THADV1C.313    
     &                  P_FIELD,extended_P_FIELD,                          THADV1C.314    
     &                  ROW_LENGTH,tot_U_ROWS,P_LEVELS,                    THADV1C.315    
     &                  EW_Halo,NS_Halo,                                   THADV1C.316    
     &                  halo_4th,halo_4th,                                 THADV1C.317    
     &                  .TRUE.)                                            THADV1C.318    
                                                                           THADV1C.319    
       ENDIF ! IF (L_SECOND)                                               THADV1C.320    
                                                                           THADV1C.321    
CL LOOP OVER P_LEVELS+1.                                                   THADV1C.322    
CL    ON 1 TO P_LEVELS PROVISIONAL VALUES OF THE FIELD ARE CALCULATED.     THADV1C.323    
CL    ON 2 TO P_LEVELS+1 THE FINAL INCREMENTS ARE CALCULATED AND ADDED     THADV1C.324    
CL    ON. THE REASON FOR THIS LOGIC IS THAT THE PROVISIONAL VALUE AT       THADV1C.325    
CL    LEVEL K+1 IS NEEDED BEFORE THE FINAL INCREMENT AT LEVEL K CAN BE     THADV1C.326    
CL    CALCULATED.                                                          THADV1C.327    
                                                                           THADV1C.328    
      DO K=1,P_LEVELS+1                                                    THADV1C.329    
CL SET TIMESTEP APPROPRIATE TO LEVEL                                       THADV1C.330    
                                                                           THADV1C.331    
        TIMESTEP = ADVECTION_TIMESTEP                                      THADV1C.332    
CL IF NOT AT P_LEVELS+1 THEN                                               THADV1C.333    
        IF(K.LE.P_LEVELS) THEN                                             THADV1C.334    
                                                                           THADV1C.335    
CL---------------------------------------------------------------------    THADV1C.336    
CL    SECTION 2.     CALCULATE COURANT NUMBER DEPENDENT NU IF IN           THADV1C.337    
CL                   FORECAST MODE.                                        THADV1C.338    
CL                   CALCULATE                                             THADV1C.339    
CL                   PROVISIONAL VALUES OF THETAL AT NEW TIME-LEVEL.       THADV1C.340    
CL---------------------------------------------------------------------    THADV1C.341    
                                                                           THADV1C.342    
C ---------------------------------------------------------------------    THADV1C.343    
CL    SECTION 2.1    SET NU TO NU_BASIC DEPENDENT ON MAX COURANT           THADV1C.344    
CL                   NUMBER.                                               THADV1C.345    
C ---------------------------------------------------------------------    THADV1C.346    
CL    IF NU_BASIC IS ZERO THEN DO NOT BOTHER TO CALCULATE NU               THADV1C.347    
          IF(.NOT.L_SECOND) THEN                                           THADV1C.348    
CL    CALCULATE COURANT NUMBER                                             THADV1C.349    
C NOTE: RS AND TRIG TERMS WILL BE INCLUDED AFTER INTERPOLATION TO P        THADV1C.350    
C       GRID.                                                              THADV1C.351    
CL    CALL UV_TO_P TO MOVE MEAN VELOCITIES ONTO P GRID                     THADV1C.352    
                                                                           THADV1C.353    
          CALL UV_TO_P(U_MEAN(START_U_REQUIRED,K),                         THADV1C.354    
     *                 NUX(START_POINT_NO_HALO,K),U_POINTS_REQUIRED,       THADV1C.355    
     *                 P_POINTS_UPDATE,ROW_LENGTH,upd_P_ROWS+1)            THADV1C.356    
                                                                           THADV1C.357    
          CALL UV_TO_P(V_MEAN(START_U_REQUIRED,K),                         THADV1C.358    
     *                 NUY(START_POINT_NO_HALO,K),U_POINTS_REQUIRED,       THADV1C.359    
     *                 P_POINTS_UPDATE,ROW_LENGTH,upd_P_ROWS+1)            THADV1C.360    
                                                                           THADV1C.361    
CL    CALCULATE NU FROM COURANT NUMBER INCLUDING TRIG AND RS TERMS.        THADV1C.362    
          DO I=START_POINT_NO_HALO,END_P_POINT_NO_HALO                     THADV1C.363    
            NUX(I,K) = NUX(I,K)*LONGITUDE_STEP_INVERSE                     THADV1C.364    
            NUY(I,K) = NUY(I,K)*LATITUDE_STEP_INVERSE                      THADV1C.365    
            SCALAR1 = TIMESTEP/(RS(I,K)*                                   THADV1C.366    
     *                RS(I,K)*(DELTA_AK(K)+DELTA_BK(K)*PSTAR_OLD(I)))      THADV1C.367    
            SCALAR2 = SEC_P_LATITUDE(I)*SCALAR1                            THADV1C.368    
            SCALAR1 = SCALAR1*SCALAR1                                      THADV1C.369    
            SCALAR2 = SCALAR2*SCALAR2                                      THADV1C.370    
            NUX(I,K) = (1. - NUX(I,K)*NUX(I,K)*SCALAR2)*NU_BASIC           THADV1C.371    
            NUY(I,K) = (1. - NUY(I,K)*NUY(I,K)*SCALAR1)*NU_BASIC           THADV1C.372    
          ENDDO                                                            THADV1C.373    
                                                                           THADV1C.374    
! Set NUX equal to minimum value along each row                            THADV1C.375    
                                                                           THADV1C.376    
          DO J=FIRST_ROW,FIRST_ROW+upd_P_ROWS-1                            THADV1C.377    
            I_start=(J-1)*ROW_LENGTH+FIRST_ROW_PT ! start and end of rpw   THADV1C.378    
            I_end=(J-1)*ROW_LENGTH+LAST_ROW_PT    ! missing out halos      THADV1C.379    
! Calculate minimum along this row                                         THADV1C.380    
*IF DEF,CRAY                                                               THADV1C.381    
            IK=ISMIN(I_end-I_start+1,NUX(I_start,K),1)                     THADV1C.382    
            SCALAR1=NUX(IK+I_start-1,K)                                    THADV1C.383    
*ELSE                                                                      THADV1C.384    
            SCALAR1=NUX(I_start,K)                                         THADV1C.385    
            DO I=I_start+1,I_end                                           THADV1C.386    
              IF (NUX(I,K) .LT. SCALAR1) SCALAR1=NUX(I,K)                  THADV1C.387    
            ENDDO                                                          THADV1C.388    
*ENDIF                                                                     THADV1C.389    
            NUX_MIN(J-FIRST_ROW+1)=SCALAR1                                 THADV1C.390    
! The indexing of NUX_MIN goes from 1..ROWS                                THADV1C.391    
          ENDDO ! J : loop over rows                                       THADV1C.392    
                                                                           THADV1C.393    
! So far we have only calculated the minimum along our local               THADV1C.394    
! part of the row. Now we must find the minimum of all the                 THADV1C.395    
! local minimums along the row                                             THADV1C.396    
          CALL GCG_RMIN(upd_P_ROWS,GC_ROW_GROUP,info,NUX_MIN)              THADV1C.397    
                                                                           THADV1C.398    
! and now set all values of NUX to the minimum along the row               THADV1C.399    
          DO J=FIRST_ROW,FIRST_ROW+upd_P_ROWS-1                            THADV1C.400    
            IF (NUX_MIN(J-FIRST_ROW+1) .LT. 0.0)                           THADV1C.401    
     &        NUX_MIN(J-FIRST_ROW+1)=0.0                                   THADV1C.402    
                                                                           THADV1C.403    
            I_start=(J-1)*ROW_LENGTH+1  ! beginning and                    THADV1C.404    
            I_end=J*ROW_LENGTH          ! end of row                       THADV1C.405    
                                                                           THADV1C.406    
            DO I=I_start,I_end                                             THADV1C.407    
              NUX(I,K)=NUX_MIN(J-FIRST_ROW+1)                              THADV1C.408    
            ENDDO                                                          THADV1C.409    
                                                                           THADV1C.410    
          ENDDO ! J : loop over rows                                       THADV1C.411    
                                                                           THADV1C.412    
! Set NUY equal to minimum value along each column                         THADV1C.413    
                                                                           THADV1C.414    
          DO J=FIRST_ROW_PT,LAST_ROW_PT                                    GPB5F403.18     
            I_start=(FIRST_ROW-1)*ROW_LENGTH+J                             THADV1C.417    
! I_start points to the beginning of column J                              THADV1C.418    
                                                                           THADV1C.419    
! Calculate the minimum along this column                                  THADV1C.420    
*IF DEF,CRAY                                                               THADV1C.421    
            IK=ISMIN(upd_P_ROWS,NUY(I_start,K),ROW_LENGTH)                 THADV1C.422    
            SCALAR1=NUY((IK-1)*ROW_LENGTH+I_start,K)                       GPB5F403.19     
*ELSE                                                                      THADV1C.424    
            I_end=I_start+(upd_P_ROWS-1)*ROW_LENGTH                        GPB5F403.20     
! I_end points to the end of column J                                      THADV1C.426    
            SCALAR1=NUY(I_start,K)                                         THADV1C.427    
            DO I=I_start+ROW_LENGTH,I_end,ROW_LENGTH                       THADV1C.428    
              IF (NUY(I,K) .LT. SCALAR1) SCALAR1=NUY(I,K)                  THADV1C.429    
            ENDDO                                                          THADV1C.430    
*ENDIF                                                                     THADV1C.431    
            NUY_MIN(J)=SCALAR1                                             THADV1C.432    
                                                                           THADV1C.433    
          ENDDO ! J : loop over columns                                    THADV1C.434    
! Once again, this is only the minimum along our local part                THADV1C.435    
! of each column. We must now find the miniumum of all the local           THADV1C.436    
! minimums along the column                                                THADV1C.437    
          CALL GCG_RMIN(ROW_LENGTH-2*EW_Halo,GC_COL_GROUP,info,            GPB5F403.21     
     &                  NUY_MIN(EW_Halo+1))                                GPB5F403.22     
                                                                           THADV1C.439    
! and now set all values of NUY to the minimum along the column            THADV1C.440    
          DO J=FIRST_ROW_PT,LAST_ROW_PT                                    GPB5F403.23     
            IF (NUY_MIN(J) .LT. 0.0) NUY_MIN(J)=0.0                        THADV1C.443    
                                                                           THADV1C.444    
            I_start=(FIRST_ROW-1)*ROW_LENGTH+J                             THADV1C.445    
            I_end=I_start+(upd_P_ROWS-1)*ROW_LENGTH                        GPB5F403.24     
                                                                           THADV1C.447    
            DO I=I_start,I_end,ROW_LENGTH                                  THADV1C.448    
              NUY(I,K)=NUY_MIN(J)                                          THADV1C.449    
            ENDDO                                                          THADV1C.450    
                                                                           THADV1C.451    
          ENDDO ! J : loop over columns                                    THADV1C.452    
                                                                           THADV1C.453    
        ENDIF  ! IF its fourth order advection                             THADV1C.454    
                                                                           THADV1C.455    
                                                                           THADV1C.456    
C ---------------------------------------------------------------------    THADV1C.457    
CL    SECTION 2.2    CALL ADV_P_GD TO OBTAIN FIRST INCREMENT DUE TO        THADV1C.458    
CL                   ADVECTION.                                            THADV1C.459    
C ---------------------------------------------------------------------    THADV1C.460    
                                                                           THADV1C.461    
CL    CALL ADV_P_GD FOR THETAL.                                            THADV1C.462    
          K1=K+1                                                           THADV1C.463    
                                                                           THADV1C.464    
          IF(K.EQ.P_LEVELS) THEN                                           THADV1C.465    
C PASS ANY THETAL VALUES AS THOSE APPARENTLY AT LEVEL K+1 AS ETADOT        THADV1C.466    
C IS SET TO ZERO BY USING ARRAY ZERO.                                      THADV1C.467    
            K1 = K-1                                                       THADV1C.468    
                                                                           THADV1C.469    
          CALL ADV_P_GD(THETAL(1,K1),THETAL(1,K),THETAL(1,K1),             THADV1C.470    
     *                  U_MEAN_COPY(1,K),V_MEAN_COPY(1,K),                 THADV1C.471    
     &                  ETADOT_MEAN(1,K),ZERO,SEC_P_LATITUDE,              THADV1C.472    
     *                  THETAL_FIRST_INC(1,K),NUX(1,K),NUY(1,K),P_FIELD,   THADV1C.473    
     *                  U_FIELD,ROW_LENGTH,                                THADV1C.474    
*CALL ARGFLDPT                                                             THADV1C.475    
     &                  TIMESTEP,LATITUDE_STEP_INVERSE,                    THADV1C.476    
     *                  LONGITUDE_STEP_INVERSE,SEC_U_LATITUDE,             THADV1C.477    
     *                  BRSP(1,K),L_SECOND,LWHITBROM,                      THADV1C.478    
     &                  T_COPY(1,K),extended_P_FIELD,extended_U_FIELD,     GSS1F403.710    
     &                  extended_address)                                  GSS1F403.711    
          ELSE IF(K.EQ.1)THEN                                              THADV1C.480    
                                                                           THADV1C.481    
C PASS ANY THETAL VALUES FOR LEVEL K-1 AS ETADOT AT LEVEL 1                THADV1C.482    
C IS SET TO ZERO BY USING ARRAY ZERO.                                      THADV1C.483    
          CALL ADV_P_GD(THETAL(1,K1),THETAL(1,K),THETAL(1,K1),             THADV1C.484    
     *                  U_MEAN_COPY(1,K),V_MEAN_COPY(1,K),ZERO,            THADV1C.485    
     *                  ETADOT_MEAN(1,K1),                                 THADV1C.486    
     *                  SEC_P_LATITUDE,THETAL_FIRST_INC(1,K),              THADV1C.487    
     *                  NUX(1,K),NUY(1,K),                                 THADV1C.488    
     *                  P_FIELD,U_FIELD,ROW_LENGTH,                        THADV1C.489    
*CALL ARGFLDPT                                                             THADV1C.490    
     &                  TIMESTEP,                                          THADV1C.491    
     *                  LATITUDE_STEP_INVERSE,LONGITUDE_STEP_INVERSE,      THADV1C.492    
     *                  SEC_U_LATITUDE,BRSP(1,K),L_SECOND,LWHITBROM,       THADV1C.493    
     &                  T_COPY(1,K),extended_P_FIELD,extended_U_FIELD,     GSS1F403.712    
     &                  extended_address)                                  GSS1F403.713    
          ELSE                                                             THADV1C.495    
          CALL ADV_P_GD(THETAL(1,K-1),THETAL(1,K),THETAL(1,K1),            THADV1C.496    
     *                  U_MEAN_COPY(1,K),V_MEAN_COPY(1,K),                 THADV1C.497    
     &                  ETADOT_MEAN(1,K),ETADOT_MEAN(1,K1),                THADV1C.498    
     *                  SEC_P_LATITUDE,THETAL_FIRST_INC(1,K),              THADV1C.499    
     *                  NUX(1,K),NUY(1,K),                                 THADV1C.500    
     *                  P_FIELD,U_FIELD,ROW_LENGTH,                        THADV1C.501    
*CALL ARGFLDPT                                                             THADV1C.502    
     &                  TIMESTEP,                                          THADV1C.503    
     *                  LATITUDE_STEP_INVERSE,LONGITUDE_STEP_INVERSE,      THADV1C.504    
     *                  SEC_U_LATITUDE,BRSP(1,K),L_SECOND,LWHITBROM,       THADV1C.505    
     &                  T_COPY(1,K),extended_P_FIELD,extended_U_FIELD,     GSS1F403.714    
     &                  extended_address)                                  GSS1F403.715    
                                                                           THADV1C.507    
          END IF                                                           THADV1C.508    
                                                                           THADV1C.509    
                                                                           THADV1C.510    
C ---------------------------------------------------------------------    THADV1C.511    
CL    SECTION 2.3    REMOVE MASS-WEIGHTING FROM INCREMENT AND ADD ONTO     THADV1C.512    
CL                   FIELD TO OBTAIN PROVISIONAL VALUE.                    THADV1C.513    
C ---------------------------------------------------------------------    THADV1C.514    
                                                                           THADV1C.515    
          DO I=1,START_POINT_NO_HALO-1                                     ADR2F402.3      
            THETAL_PROV(I,K) = 0.0                                         ADR2F402.4      
          ENDDO                                                            ADR2F402.5      
                                                                           ADR2F402.6      
          DO I=START_POINT_NO_HALO,END_P_POINT_NO_HALO                     THADV1C.516    
            SCALAR1 = RS(I,K)*RS(I,K)                                      THADV1C.517    
     *                      *(DELTA_AK(K)+DELTA_BK(K)*PSTAR_OLD(I))        THADV1C.518    
            THETAL_FIRST_INC(I,K)=THETAL_FIRST_INC(I,K)/SCALAR1            THADV1C.519    
            THETAL_PROV(I,K) = THETAL(I,K)- THETAL_FIRST_INC(I,K)          THADV1C.520    
          ENDDO                                                            THADV1C.521    
                                                                           THADV1C.522    
          DO I=END_P_POINT_NO_HALO+1,P_FIELD                               ADR2F402.7      
            THETAL_PROV(I,K) = 0.0                                         ADR2F402.8      
          ENDDO                                                            ADR2F402.9      
                                                                           ADR2F402.10     
*IF DEF,GLOBAL                                                             THADV1C.523    
CL   GLOBAL MODEL CALCULATE PROVISIONAL POLAR VALUE.                       THADV1C.524    
          IF (at_top_of_LPG) THEN                                          THADV1C.525    
! North Pole                                                               THADV1C.526    
            DO I=TOP_ROW_START,TOP_ROW_START+ROW_LENGTH-1                  THADV1C.527    
              THETAL_PROV(I,K) = THETAL(I,K)                               THADV1C.528    
              THETAL_FIRST_INC(I,K) = -THETAL_FIRST_INC(I,K)/              THADV1C.529    
     &                         (RS(I,K)*RS(I,K)*                           THADV1C.530    
     &                         (DELTA_AK(K)+DELTA_BK(K)*PSTAR_OLD(I)))     THADV1C.531    
            ENDDO                                                          THADV1C.532    
          ENDIF                                                            THADV1C.533    
                                                                           THADV1C.534    
          IF (at_base_of_LPG) THEN                                         THADV1C.535    
! South Pole                                                               THADV1C.536    
            DO I=P_BOT_ROW_START,P_BOT_ROW_START+ROW_LENGTH-1              THADV1C.537    
              THETAL_PROV(I,K) = THETAL(I,K)                               THADV1C.538    
              THETAL_FIRST_INC(I,K) = -THETAL_FIRST_INC(I,K)/              THADV1C.539    
     &                         (RS(I,K)*RS(I,K)*                           THADV1C.540    
     &                         (DELTA_AK(K)+DELTA_BK(K)*PSTAR_OLD(I)))     THADV1C.541    
            ENDDO                                                          THADV1C.542    
          ENDIF                                                            THADV1C.543    
                                                                           THADV1C.544    
*ELSE                                                                      THADV1C.545    
CL    LIMITED AREA MODEL THEN SET PROVISIONAL VALUES ON BOUNDARIES         THADV1C.546    
CL    TO FIELD VALUES AT OLD TIME LEVEL.                                   THADV1C.547    
      IF (at_top_of_LPG) THEN                                              THADV1C.548    
        DO I=TOP_ROW_START,TOP_ROW_START+ROW_LENGTH-1                      THADV1C.549    
          THETAL_PROV(I,K) = THETAL(I,K)                                   THADV1C.550    
        ENDDO                                                              THADV1C.551    
      ENDIF                                                                THADV1C.552    
      IF (at_base_of_LPG) THEN                                             THADV1C.553    
        DO I=P_BOT_ROW_START,P_BOT_ROW_START+ROW_LENGTH-1                  THADV1C.554    
          THETAL_PROV(I,K) = THETAL(I,K)                                   THADV1C.555    
        ENDDO                                                              THADV1C.556    
      ENDIF                                                                THADV1C.557    
*ENDIF                                                                     THADV1C.558    
        END IF                                                             THADV1C.559    
CL END CONDITIONAL ON LEVEL BEING LESS THAN P_LEVELS+1                     THADV1C.560    
      ENDDO                                                                THADV1C.561    
                                                                           THADV1C.562    
*IF DEF,GLOBAL                                                             THADV1C.563    
                                                                           THADV1C.564    
      CALL POLAR(THETAL_PROV,THETAL_FIRST_INC,THETAL_FIRST_INC,            THADV1C.565    
*CALL ARGFLDPT                                                             THADV1C.566    
     &           P_FIELD,P_FIELD,P_FIELD,                                  THADV1C.567    
     &           TOP_ROW_START,P_BOT_ROW_START,                            THADV1C.568    
     &           ROW_LENGTH,P_LEVELS)                                      THADV1C.569    
*ENDIF                                                                     THADV1C.570    
                                                                           THADV1C.571    
      IF (L_SECOND) THEN                                                   THADV1C.572    
                                                                           THADV1C.573    
! Do a halo update on the THETAL_PROV array                                THADV1C.574    
! that has just been calculated                                            THADV1C.575    
        CALL SWAPBOUNDS(THETAL_PROV,ROW_LENGTH,tot_P_ROWS,                 THADV1C.576    
     &                  EW_Halo,NS_Halo,P_LEVELS)                          THADV1C.577    
!        CALL SET_SIDES(THETAL_PROV,P_FIELD,ROW_LENGTH,P_LEVELS,           THADV1C.578    
!     &                 fld_type_p)                                        THADV1C.579    
                                                                           THADV1C.580    
      ELSE  ! fourth order advection                                       THADV1C.581    
                                                                           THADV1C.582    
! Copy THETAL_PROV into T_COPY which has double halos for fourth           THADV1C.583    
! order advection, and do swap to fill these halos                         THADV1C.584    
        CALL COPY_FIELD(THETAL_PROV,T_COPY,                                THADV1C.585    
     &                  P_FIELD,extended_P_FIELD,                          THADV1C.586    
     &                  ROW_LENGTH,tot_P_ROWS,P_LEVELS,                    THADV1C.587    
     &                  EW_Halo,NS_Halo,                                   THADV1C.588    
     &                  halo_4th,halo_4th,                                 THADV1C.589    
     &                  .TRUE.)                                            THADV1C.590    
                                                                           THADV1C.591    
      ENDIF                                                                THADV1C.592    
                                                                           THADV1C.593    
! Set up OMEGA_P array                                                     THADV1C.594    
! (Was SECTION 3.2):                                                       THADV1C.595    
!    SECTION 3.2    INTERPOLATE OMEGA TO P GRID AND CALCULATE              THADV1C.596    
!                   REMAINING TERM IN ADVECTION EQUATION.                  THADV1C.597    
!                   CALCULATE TOTAL MASS-WEIGHTED INCREMENT TO FIELD.      THADV1C.598    
                                                                           THADV1C.599    
            DO K1=1,P_LEVELS                                               THADV1C.600    
              CALL UV_TO_P(OMEGA(START_U_REQUIRED,K1),                     THADV1C.601    
     &                     OMEGA_P(START_POINT_NO_HALO,K1),                THADV1C.602    
     &                     U_POINTS_REQUIRED,                              THADV1C.603    
     &                     P_POINTS_UPDATE,ROW_LENGTH,upd_P_ROWS+1)        THADV1C.604    
*IF DEF,GLOBAL                                                             THADV1C.605    
              IF (at_top_of_LPG) THEN                                      THADV1C.606    
                DO I=TOP_ROW_START,TOP_ROW_START+ROW_LENGTH-1              THADV1C.607    
                  OMEGA_P(I,K1)=0.                                         THADV1C.608    
                ENDDO                                                      THADV1C.609    
              ENDIF                                                        THADV1C.610    
                                                                           THADV1C.611    
              IF (at_base_of_LPG) THEN                                     THADV1C.612    
                DO I=P_BOT_ROW_START,P_BOT_ROW_START+ROW_LENGTH-1          THADV1C.613    
                  OMEGA_P(I,K1)=0.                                         THADV1C.614    
                ENDDO                                                      THADV1C.615    
              ENDIF                                                        THADV1C.616    
*ENDIF                                                                     THADV1C.617    
            ENDDO                                                          THADV1C.618    
*IF DEF,GLOBAL                                                             THADV1C.619    
                                                                           THADV1C.620    
              CALL POLAR(OMEGA_P,OMEGA_P,OMEGA_P,                          THADV1C.621    
*CALL ARGFLDPT                                                             THADV1C.622    
     &                   P_FIELD,P_FIELD,P_FIELD,                          THADV1C.623    
     &                   START_POINT_NO_HALO,                              THADV1C.624    
     &                   END_P_POINT_NO_HALO-ROW_LENGTH+1,                 THADV1C.625    
     &                   ROW_LENGTH,P_LEVELS)                              THADV1C.626    
*ENDIF                                                                     THADV1C.627    
                                                                           THADV1C.628    
CL BEGIN CONDITIONAL ON LEVEL BEING GREATER THAN 1                         THADV1C.629    
                                                                           THADV1C.630    
      DO K=1,P_LEVELS+1                                                    THADV1C.631    
        IF(K.GT.1) THEN                                                    THADV1C.632    
CL---------------------------------------------------------------------    THADV1C.633    
CL    SECTION 3.     ALL WORK IN THIS SECTION PERFORMED AT LEVEL-1.        THADV1C.634    
CL                   CALCULATE SECOND INCREMENT DUE TO ADVECTION.          THADV1C.635    
CL                   CALCULATE TOTAL INCREMENT TO FIELD AND FILTER         THADV1C.636    
CL                   WHERE NECESSARY THEN UPDATE FIELD.                    THADV1C.637    
CL                   THE POLAR INCREMENTS ARE THEN CALCULATED AND ADDED    THADV1C.638    
CL                   ON BY CALLING POLAR.                                  THADV1C.639    
CL---------------------------------------------------------------------    THADV1C.640    
                                                                           THADV1C.641    
         TIMESTEP = ADVECTION_TIMESTEP                                     THADV1C.642    
                                                                           THADV1C.643    
         CONST1 = R/(CP*CP)*TIMESTEP                                       THADV1C.644    
C ---------------------------------------------------------------------    THADV1C.645    
CL    SECTION 3.1    CALL ADV_P_GD TO OBTAIN SECOND INCREMENT DUE TO       THADV1C.646    
CL                   ADVECTION.                                            THADV1C.647    
C ---------------------------------------------------------------------    THADV1C.648    
                                                                           THADV1C.649    
CL    CALL ADV_P_GD FOR THETAL.                                            THADV1C.650    
          K1=K-1                                                           THADV1C.651    
C K1 HOLDS K-1.                                                            THADV1C.652    
          IF(K.GT.P_LEVELS) THEN                                           THADV1C.653    
C THE ZERO VERTICAL FLUX AT THE TOP IS ENSURED BY PASSING ETADOT AS        THADV1C.654    
C ZERO.                                                                    THADV1C.655    
                                                                           THADV1C.656    
          CALL ADV_P_GD(THETAL_PROV(1,K-2),THETAL_PROV(1,K-1),             THADV1C.657    
     *                  THETAL_PROV(1,K-2),                                THADV1C.658    
     *                  U_MEAN_COPY(1,K1),V_MEAN_COPY(1,K1),               THADV1C.659    
     &                  ETADOT_MEAN(1,K-1),ZERO,SEC_P_LATITUDE,            THADV1C.660    
     *                  THETAL_SECOND_INC,NUX(1,K-1),NUY(1,K-1),P_FIELD,   THADV1C.661    
     *                  U_FIELD,ROW_LENGTH,                                THADV1C.662    
*CALL ARGFLDPT                                                             THADV1C.663    
     &                  TIMESTEP,LATITUDE_STEP_INVERSE,                    THADV1C.664    
     *                  LONGITUDE_STEP_INVERSE,SEC_U_LATITUDE,             THADV1C.665    
     *                  BRSP(1,K-1),L_SECOND,LWHITBROM,                    THADV1C.666    
     &                  T_COPY(1,K-1),extended_P_FIELD,extended_U_FIELD,   GSS1F403.716    
     &                  extended_address)                                  GSS1F403.717    
                                                                           THADV1C.668    
          ELSE IF(K.EQ.2) THEN                                             THADV1C.669    
C THE ZERO VERTICAL FLUX AT THE BOTTOM IS ENSURED BY PASSING ETADOT AS     THADV1C.670    
C ZERO.                                                                    THADV1C.671    
          CALL ADV_P_GD(THETAL_PROV(1,K),THETAL_PROV(1,K-1),               THADV1C.672    
     *                  THETAL_PROV(1,K),                                  THADV1C.673    
     *                  U_MEAN_COPY(1,K1),V_MEAN_COPY(1,K1),ZERO,          THADV1C.674    
     *                  ETADOT_MEAN(1,K),                                  THADV1C.675    
     *                  SEC_P_LATITUDE,THETAL_SECOND_INC,                  THADV1C.676    
     *                  NUX(1,K-1),NUY(1,K-1),                             THADV1C.677    
     *                  P_FIELD,U_FIELD,ROW_LENGTH,                        THADV1C.678    
*CALL ARGFLDPT                                                             THADV1C.679    
     &                  TIMESTEP,                                          THADV1C.680    
     *                  LATITUDE_STEP_INVERSE,LONGITUDE_STEP_INVERSE,      THADV1C.681    
     *                  SEC_U_LATITUDE,                                    THADV1C.682    
     *                  BRSP(1,K-1),L_SECOND,LWHITBROM,                    THADV1C.683    
     &                  T_COPY(1,K-1),extended_P_FIELD,extended_U_FIELD,   GSS1F403.718    
     &                  extended_address)                                  GSS1F403.719    
          ELSE                                                             THADV1C.685    
                                                                           THADV1C.686    
          CALL ADV_P_GD(THETAL_PROV(1,K-2),THETAL_PROV(1,K-1),             THADV1C.687    
     *                  THETAL_PROV(1,K),                                  THADV1C.688    
     *                  U_MEAN_COPY(1,K1),V_MEAN_COPY(1,K1),               THADV1C.689    
     &                  ETADOT_MEAN(1,K-1),ETADOT_MEAN(1,K),               THADV1C.690    
     *                  SEC_P_LATITUDE,THETAL_SECOND_INC,                  THADV1C.691    
     *                  NUX(1,K-1),NUY(1,K-1),                             THADV1C.692    
     *                  P_FIELD,U_FIELD,ROW_LENGTH,                        THADV1C.693    
*CALL ARGFLDPT                                                             THADV1C.694    
     &                  TIMESTEP,                                          THADV1C.695    
     *                  LATITUDE_STEP_INVERSE,LONGITUDE_STEP_INVERSE,      THADV1C.696    
     *                  SEC_U_LATITUDE,                                    THADV1C.697    
     *                  BRSP(1,K-1),L_SECOND,LWHITBROM,                    THADV1C.698    
     &                  T_COPY(1,K-1),extended_P_FIELD,extended_U_FIELD,   GSS1F403.720    
     &                  extended_address)                                  GSS1F403.721    
          END IF                                                           THADV1C.700    
                                                                           THADV1C.701    
                                                                           THADV1C.702    
C TOTAL MASS-WEIGHTED HORIZONTAL AND VERTICAL INCREMENTS ARE CALCULATED    THADV1C.703    
C SEPARATELY.                                                              THADV1C.704    
                                                                           THADV1C.705    
          IF(K.LT.Q_LEVELS+2) THEN                                         THADV1C.706    
            DO I=START_POINT_NO_HALO,END_P_POINT_NO_HALO                   THADV1C.707    
                                                                           THADV1C.708    
              PK  = AKH(K)  + BKH(K) *PSTAR(I)                             THADV1C.709    
              PK1 = AKH(K1) + BKH(K1)*PSTAR(I)  !  K1 = K-1                THADV1C.710    
              P_EXNER_FULL = P_EXNER_C                                     THADV1C.711    
     *        (P_EXNER(I,K),P_EXNER(I,K1),PK,PK1,KAPPA)                    THADV1C.712    
                                                                           THADV1C.713    
              THETAL_INCREMENT(I,K1) = .5*(THETAL_SECOND_INC(I) +          THADV1C.714    
     *                       THETAL_FIRST_INC(I,K-1)*RS(I,K1)*RS(I,K1)     THADV1C.715    
     *                      *(DELTA_AK(K1)+DELTA_BK(K1)*PSTAR(I)))         THADV1C.716    
     *                      -(LC*QCL(I,K1)+LC_LF*QCF(I,K1))*CONST1*        THADV1C.717    
     &                       OMEGA_P(I,K1)/((AK(K1)+BK(K1)*PSTAR(I))       THADV1C.718    
     *                       *(P_EXNER_FULL))                              THADV1C.719    
                                                                           THADV1C.720    
            ENDDO                                                          THADV1C.721    
          ELSE                                                             THADV1C.722    
            DO I=START_POINT_NO_HALO,END_P_POINT_NO_HALO                   THADV1C.723    
              THETAL_INCREMENT(I,K1) = .5*(THETAL_SECOND_INC(I) +          THADV1C.724    
     *                       THETAL_FIRST_INC(I,K-1)*RS(I,K1)*RS(I,K1)     THADV1C.725    
     *                      *(DELTA_AK(K1)+DELTA_BK(K1)*PSTAR(I)))         THADV1C.726    
            ENDDO                                                          THADV1C.727    
          END IF                                                           THADV1C.728    
                                                                           THADV1C.729    
C ---------------------------------------------------------------------    THADV1C.730    
CL    SECTION 3.3    IF GLOBAL MODEL CALCULATE POLAR INCREMENTS.           THADV1C.731    
CL                   IF LIMITED AREA MASS-WEIGHT BOUNDARIES.               THADV1C.732    
C ---------------------------------------------------------------------    THADV1C.733    
                                                                           THADV1C.734    
CL    GLOBAL MODEL CALCULATE POLAR INCREMENT.                              THADV1C.735    
CL    CALCULATE MERIDIONAL FLUX AROUND POLES BY ADDING THE TWO             THADV1C.736    
CL    INCREMENTS AND ALSO MASS-WEIGHTING POLAR FIELDS.                     THADV1C.737    
C NEGATIVE SIGN BEFORE FIRST INCS IS DUE TO THEIR SIGN HAVING BEEN         THADV1C.738    
C CHANGED PRIOR TO THE CALCULATION OF THE INTERMEDIATE VALUE.              THADV1C.739    
          IF (at_top_of_LPG) THEN                                          THADV1C.740    
! Northern boundary/pole                                                   THADV1C.741    
            IF (K.LT.Q_LEVELS+2) THEN                                      THADV1C.742    
              DO I=TOP_ROW_START,TOP_ROW_START+ROW_LENGTH-1                THADV1C.743    
                SCALAR1 = RS(I,K1)*RS(I,K1)*                               THADV1C.744    
     &            (DELTA_AK(K1)+DELTA_BK(K1)*PSTAR(I))                     THADV1C.745    
*IF DEF,GLOBAL                                                             THADV1C.746    
                PK  = AKH(K)  + BKH(K) *PSTAR(I)                           THADV1C.747    
                PK1 = AKH(K1) + BKH(K1)*PSTAR(I)  !  K1 = K-1              THADV1C.748    
                P_EXNER_FULL = P_EXNER_C(P_EXNER(I,K),                     THADV1C.749    
     &                                   P_EXNER(I,K1),PK,PK1,KAPPA)       THADV1C.750    
                                                                           THADV1C.751    
                THETAL_INCREMENT(I,K1) = -.5*(THETAL_SECOND_INC(I)         THADV1C.752    
     &                           - THETAL_FIRST_INC(I,K-1)*SCALAR1)        THADV1C.753    
     &                           +(LC*QCL(I,K1)+LC_LF*QCF(I,K1))*CONST1*   THADV1C.754    
     &                           OMEGA_P(I,K1)/((AK(K1)+BK(K1)*PSTAR(I))   THADV1C.755    
     &                           *P_EXNER_FULL)                            THADV1C.756    
*ENDIF                                                                     THADV1C.757    
                THETAL(I,K1) = THETAL(I,K1)*SCALAR1                        THADV1C.758    
              ENDDO                                                        THADV1C.759    
            ELSE  ! (IF K.GE.Q_LEVELS+2)                                   THADV1C.760    
              DO I=TOP_ROW_START,TOP_ROW_START+ROW_LENGTH-1                THADV1C.761    
                SCALAR1 = RS(I,K1)*RS(I,K1)*                               THADV1C.762    
     &                    (DELTA_AK(K1)+DELTA_BK(K1)*PSTAR(I))             THADV1C.763    
*IF DEF,GLOBAL                                                             THADV1C.764    
                THETAL_INCREMENT(I,K1) = -.5*(THETAL_SECOND_INC(I)         THADV1C.765    
     &                              - THETAL_FIRST_INC(I,K-1)*SCALAR1)     THADV1C.766    
*ENDIF                                                                     THADV1C.767    
                THETAL(I,K1) = THETAL(I,K1)*SCALAR1                        THADV1C.768    
              ENDDO                                                        THADV1C.769    
            ENDIF ! (K.LT.Q_LEVELS+2)                                      THADV1C.770    
          ENDIF ! (attop)                                                  THADV1C.771    
                                                                           THADV1C.772    
          IF (at_base_of_LPG) THEN                                         THADV1C.773    
! Southern boundary/pole                                                   THADV1C.774    
            IF (K.LT.Q_LEVELS+2) THEN                                      THADV1C.775    
              DO I=P_BOT_ROW_START,P_BOT_ROW_START+ROW_LENGTH-1            THADV1C.776    
                SCALAR2 = RS(I,K1)*RS(I,K1)*                               THADV1C.777    
     &            (DELTA_AK(K1)+DELTA_BK(K1)*PSTAR(I))                     THADV1C.778    
*IF DEF,GLOBAL                                                             THADV1C.779    
                PK  = AKH(K)  + BKH(K) *PSTAR(I)                           THADV1C.780    
                PK1 = AKH(K1) + BKH(K1)*PSTAR(I)  !  K1 = K-1              THADV1C.781    
                P_EXNER_FULL = P_EXNER_C(P_EXNER(I,K),                     THADV1C.782    
     &                                   P_EXNER(I,K1),PK,PK1,KAPPA)       THADV1C.783    
                                                                           THADV1C.784    
                THETAL_INCREMENT(I,K1) = -.5*(THETAL_SECOND_INC(I)         THADV1C.785    
     &                           - THETAL_FIRST_INC(I,K-1)*SCALAR2)        THADV1C.786    
     &                           +(LC*QCL(I,K1)+LC_LF*QCF(I,K1))*CONST1*   THADV1C.787    
     &                           OMEGA_P(I,K1)/((AK(K1)+BK(K1)*PSTAR(I))   THADV1C.788    
     &                           *P_EXNER_FULL)                            THADV1C.789    
*ENDIF                                                                     THADV1C.790    
                THETAL(I,K1) = THETAL(I,K1)*SCALAR2                        THADV1C.791    
              ENDDO                                                        THADV1C.792    
            ELSE  ! (IF K.GE.Q_LEVELS+2)                                   THADV1C.793    
              DO I=P_BOT_ROW_START,P_BOT_ROW_START+ROW_LENGTH-1            THADV1C.794    
                SCALAR2 = RS(I,K1)*RS(I,K1)*                               THADV1C.795    
     &                    (DELTA_AK(K1)+DELTA_BK(K1)*PSTAR(I))             THADV1C.796    
                THETAL(I,K1) = THETAL(I,K1)*SCALAR2                        THADV1C.797    
*IF DEF,GLOBAL                                                             THADV1C.798    
                THETAL_INCREMENT(I,K1) = -.5*(THETAL_SECOND_INC(I)         THADV1C.799    
     &                              - THETAL_FIRST_INC(I,K-1)*SCALAR2)     THADV1C.800    
*ENDIF                                                                     THADV1C.801    
              ENDDO                                                        THADV1C.802    
            ENDIF ! (K.LT.Q_LEVELS+2)                                      THADV1C.803    
          ENDIF ! (atbase)                                                 THADV1C.804    
                                                                           THADV1C.805    
CL END CONDITIONAL LEVEL GREATER THAN ONE                                  THADV1C.806    
        END IF                                                             THADV1C.807    
                                                                           THADV1C.808    
CL END LOOP OVER P_LEVELS+1                                                THADV1C.809    
      enddo                                                                THADV1C.810    
                                                                           THADV1C.811    
CL---------------------------------------------------------------------    THADV1C.812    
CL    SECTION 4      IF GLOBAL MODEL THEN FILTER INCREMENTS AND            THADV1C.813    
CL                   UPDATE POLAR VALUES BY CALLING POLAR.                 THADV1C.814    
CL                   UPDATE ALL OTHER VALUES.                              THADV1C.815    
CL---------------------------------------------------------------------    THADV1C.816    
                                                                           THADV1C.817    
*IF DEF,GLOBAL                                                             THADV1C.818    
                                                                           THADV1C.819    
C ---------------------------------------------------------------------    THADV1C.820    
CL    SECTION 4.1    CALL FILTER TO DO FILTERING.                          THADV1C.821    
C ---------------------------------------------------------------------    THADV1C.822    
                                                                           THADV1C.823    
C SET FILTER_SPACE WHICH IS ROW_LENGTH+2 TIMES THE NUMBER OF ROWS TO       THADV1C.824    
C BE FILTERED.                                                             THADV1C.825    
                                                                           THADV1C.826    
      FILTER_SPACE = (ROW_LENGTH+2)*(NORTHERN_FILTERED_P_ROW-1+            THADV1C.827    
     *                tot_P_ROWS-SOUTHERN_FILTERED_P_ROW)                  THADV1C.828    
CL    CALL FILTER FOR THETAL INCREMENTS                                    THADV1C.829    
                                                                           THADV1C.830    
      CALL FILTER(THETAL_INCREMENT,P_FIELD,P_LEVELS,                       APB0F402.47     
     &            FILTER_SPACE,ROW_LENGTH,                                 APB0F402.48     
*CALL ARGFLDPT                                                             APB0F402.49     
     &            FILTER_WAVE_NUMBER_P_ROWS,TRIGS,IFAX,                    APB0F402.50     
     *            NORTHERN_FILTERED_P_ROW,SOUTHERN_FILTERED_P_ROW)         THADV1C.834    
                                                                           THADV1C.835    
C ---------------------------------------------------------------------    THADV1C.836    
CL    SECTION 4.2    CALL POLAR TO UPDATE POLAR VALUES                     THADV1C.837    
C ---------------------------------------------------------------------    THADV1C.838    
                                                                           THADV1C.839    
      CALL POLAR(THETAL,THETAL_INCREMENT,THETAL_INCREMENT,                 THADV1C.840    
*CALL ARGFLDPT                                                             THADV1C.841    
     &           P_FIELD,P_FIELD,P_FIELD,                                  THADV1C.842    
     &           TOP_ROW_START,P_BOT_ROW_START,                            THADV1C.843    
     &           ROW_LENGTH,P_LEVELS)                                      THADV1C.844    
                                                                           THADV1C.845    
*ENDIF                                                                     THADV1C.846    
C ---------------------------------------------------------------------    THADV1C.847    
CL    SECTION 4.3    UPDATE ALL OTHER POINTS.                              THADV1C.848    
C   OUTPUT IS MASS-WEIGHTED.                                               THADV1C.849    
C   INCREMENTS ARE ALREADY MASS-WEIGHTED                                   THADV1C.850    
C ---------------------------------------------------------------------    THADV1C.851    
                                                                           THADV1C.852    
      DO K=1,P_LEVELS                                                      THADV1C.853    
C UPDATE THETAL.                                                           THADV1C.854    
CFPP$ SELECT(CONCUR)                                                       THADV1C.855    
        DO I= START_POINT_NO_HALO,END_P_POINT_NO_HALO                      THADV1C.856    
          THETAL(I,K)=THETAL(I,K)*RS(I,K)*RS(I,K)*                         THADV1C.857    
     &        (DELTA_AK(K)+DELTA_BK(K)*PSTAR(I))-THETAL_INCREMENT(I,K)     THADV1C.858    
        ENDDO                                                              THADV1C.859    
      ENDDO                                                                THADV1C.860    
                                                                           THADV1C.861    
CL    END OF ROUTINE TH_ADV                                                THADV1C.862    
                                                                           THADV1C.863    
      RETURN                                                               THADV1C.864    
      END                                                                  THADV1C.865    
*ENDIF                                                                     THADV1C.866