1C Copyright 1981-2016 ECMWF. 2C 3C This software is licensed under the terms of the Apache Licence 4C Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 5C 6C In applying this licence, ECMWF does not waive the privileges and immunities 7C granted to it by virtue of its status as an intergovernmental organisation 8C nor does it submit to any jurisdiction. 9C 10 11 INTEGER FUNCTION IGGLAT(KLAT, PGAUSS, KPR, KERR) 12C 13C----> 14C**** *IGGLAT* 15C 16C PURPOSE 17C _______ 18C 19C Compute Gaussian lines of latitude for a given truncation. 20C 21C 22C INTERFACE 23C _________ 24C 25C IERR = IGGLAT(KLAT, PGAUSS, KPR, KERR) 26C 27C 28C Input parameters 29C ________________ 30C 31C KLAT - The number of lines of latitude from pole to pole. 32C 33C KPR - Debug print switch: 34C 0, No debugging output. 35C 1, Produce debugging output. 36C 37C KERR - Error control flag. 38C -ve, No error message. Return error code. 39C 0 , Hard failure with error message. 40C +ve, Print error message. Return error code. 41C 42C 43C Output parameters 44C ________________ 45C 46C PGAUSS - The Gaussian lines of latitude for truncation KLAT/2 47C 48C 49C Return value 50C ____________ 51C 52C 0 - OK 53C 24701 - The calculation of Gaussian lines of latitude failed. 54C 55C 56C Common block usage 57C __________________ 58C 59C None 60C 61C 62C EXTERNALS 63C _________ 64C 65C INTLOG(R) - Logs messages. 66C IGBESS - This routine is used to generate initial 67C approximations to the Gaussian latitudes. 68C 69C 70C METHOD 71C ______ 72C 73C IGBESS is used to provide zeros of the Bessel function J0, 74C which are used as starting approximations to the Gaussian 75C latitudes. Newton iteration is used to generate the latitudes 76C from these approximations. 77C 78C Calculated latitudes are stored internally and reused if a 79C later call asks for the same calculation. 80C 81C 82C REFERENCE 83C _________ 84C 85C None 86C 87C 88C COMMENTS 89C ________ 90C 91C This routine is adapted from that in the old Marsint library. 92C The interface and the variable names have been modified. 93C 94C 95C AUTHOR 96C ______ 97C 98C K. Fielding *ECMWF* Oct 1993 99C 100C 101C MODIFICATIONS 102C _____________ 103C 104C J.Chambers ECMWF Aug 1998 105C Add reuse of calculated latitudes. 106C 107C----< 108C ------------------------------------------------------- 109C* Section 0. Definition of variables. 110C ------------------------------------------------------- 111C 112 IMPLICIT NONE 113C 114#include "jparams.h" 115#include "parim.h" 116C 117C Function arguments 118C 119 INTEGER KLAT, KPR, KERR 120 REAL PGAUSS (*) 121C 122C Parameters 123C 124 INTEGER JPROUTINE 125 PARAMETER (JPROUTINE = 24700) 126C ------------------------------------------------------- 127C The convergence criteria is machine dependent 128C ------------------------------------------------------- 129C 130 REAL PACRCY 131#ifdef REAL_8 132 PARAMETER (PACRCY = 1.0E-14) 133#else 134 PARAMETER (PACRCY = 1.0E-7) 135#endif 136C 137C Local variables 138C 139 LOGICAL LDEBUG 140 INTEGER ITRUNC, IERR 141 INTEGER JLAT, JITER, JLN, LOOP, JOLDLAT 142 REAL ZRADDEG, ZCON, ZLAT, ZROOT 143 REAL ZKM1, ZKM2, ZLN, ZFUNC, ZDERIV, ZMOVE 144 REAL ROLDGS 145 DIMENSION ROLDGS(JPGTRUNC*2) 146C 147 SAVE JOLDLAT, ROLDGS 148C 149C Externals 150C 151 INTEGER IGBESS 152C 153C ------------------------------------------------------- 154C* Section 1. Set constants and get initial approximation. 155C ------------------------------------------------------- 156C 157 100 CONTINUE 158C 159 ZFUNC = 0. 160 IGGLAT = 0 161 LDEBUG = ( KPR.GE.1 ) 162C 163 IF( LDEBUG ) THEN 164 CALL INTLOG(JP_DEBUG,'IGGLAT: Section 1.',JPQUIET) 165 CALL INTLOG(JP_DEBUG, 166 X 'IGGLAT: No.lines lat from pole to pole = ',KLAT) 167 ENDIF 168C 169C Using previously calculated values if truncation is the same. 170C 171 IF( KLAT.EQ.JOLDLAT ) THEN 172 IF( LDEBUG ) CALL INTLOG(JP_DEBUG, 173 X 'IGGLAT: Using previously calculated values',JPQUIET) 174 DO LOOP = 1, KLAT 175 PGAUSS(LOOP) = ROLDGS(LOOP) 176 ENDDO 177 GOTO 900 178 ENDIF 179C 180 JOLDLAT = KLAT 181C 182 ZRADDEG = 180.0/PPI 183C 184 ZCON = (PPONE - (PPTWO / PPI) ** 2) * PPQUART 185C 186 ZLAT = KLAT 187 ITRUNC = KLAT/2 188C 189 IERR = IGBESS(ITRUNC, PGAUSS, KPR, KERR) 190C 191 IF( IERR.GT.0 ) THEN 192 IGGLAT = IERR 193 GOTO 900 194 ENDIF 195C 196C ------------------------------------------------------- 197C* Section 2. Compute abscissae 198C ------------------------------------------------------- 199C 200 200 CONTINUE 201C 202 IF( LDEBUG ) CALL INTLOG(JP_DEBUG,'IGGLAT: Section 2.',JPQUIET) 203C 204 DO 240 JLAT = 1, ITRUNC 205C 206C First approximation for ZROOT 207C 208 ZROOT = COS(PGAUSS(JLAT) / SQRT( (ZLAT+PPHALF)**2 + ZCON) ) 209C 210C Perform loop of Newton iterations 211C 212 DO 220 JITER = 1, JPMAXITER 213C 214 ZKM2 = PPONE 215 ZKM1 = ZROOT 216C 217C Compute Legendre polynomial 218C 219 DO 210 JLN = 2, KLAT 220C 221 ZLN = JLN 222C 223 ZFUNC = ( (PPTWO * ZLN - PPONE) * ZROOT * ZKM1 - 224 1 (ZLN - PPONE) * ZKM2) / ZLN 225C 226 ZKM2 = ZKM1 227 ZKM1 = ZFUNC 228C 229 210 CONTINUE 230C 231C Perform Newton iteration 232C 233 ZDERIV = (ZLAT * (ZKM2 - ZROOT * ZFUNC) ) / 234 1 (PPONE - ZROOT ** 2) 235C 236 ZMOVE = ZFUNC / ZDERIV 237 ZROOT = ZROOT - ZMOVE 238C 239C Leave iteration loop when sufficient accuracy achieved. 240C 241 IF( ABS(ZMOVE).LE.PACRCY ) GOTO 230 242C 243 220 CONTINUE 244C 245C Routine fails if no convergence after JPMAXITER iterations. 246C 247 IGGLAT = JPROUTINE + 1 248C 249 IF( KERR.GE.0 ) CALL INTLOG(JP_ERROR, 250 X 'IGGLAT: Calculation of Gaussian lats failed.',JPQUIET) 251C 252 IF( KERR.EQ.0 ) CALL INTLOG(JP_FATAL, 253 X 'IGGLAT: interpolation failed.',IGGLAT) 254C 255 GOTO 900 256C 257 230 CONTINUE 258C 259C* Set North and South values using symmetry. 260C 261 PGAUSS(JLAT) = ASIN(ZROOT) * ZRADDEG 262 PGAUSS(KLAT + 1 - JLAT) = - PGAUSS(JLAT) 263C 264 240 CONTINUE 265C 266 IF( KLAT.NE.(ITRUNC*2) ) PGAUSS(ITRUNC + 1) = PPZERO 267C 268C Store calculated values for re-use 269C 270 DO LOOP = 1, KLAT 271 ROLDGS(LOOP) = PGAUSS(LOOP) 272 ENDDO 273C 274C ------------------------------------------------------- 275C* Section 9. Return to calling routine. 276C ------------------------------------------------------- 277C 278 900 CONTINUE 279C 280 IF( LDEBUG ) CALL INTLOG(JP_DEBUG,'IGGLAT: Section 9.',JPQUIET) 281C 282 RETURN 283 END 284