SUBROUTINE WACCMC2_ATMSTATE_INTERP(NTIMES,NLATITUDES,NLEVELS, &Z_WACCM,P_WACCM,ZGEO_WACCM,T_WACCM, &H2O_WACCM,CO2_WACCM,O3_WACCM,N2O_WACCM,CH4_WACCM,O_WACCM, &O2_WACCM,NO_WACCM,NO2_WACCM,OH_WACCM,N2_WACCM,HO2_WACCM, &CFCL3_WACCM,CF2CL2_WACCM,H_WACCM,O1D_WACCM, &MR,MGM,IDG,RE,XLAT,ZR,PR,TR,Q,CO2,CO2_FIX) ************************************************************************ * WACCM_ATMSTATE CREATES ATMOSPHERIC STATE VECTORS FROM THE WACCM * MODEL DATABASE ARRAYS. THERE IS ONE FILE PER MONTH THROUGHOUT THE * SABER MISSION. * * ZGEO_WACCM IS GEOPOTENTIAL ALTITUDE IN METERS WE MUST CONVERT THAT * TO GEOMETRIC ALTITUDE IN KM. WE TAKE ALL PROFILES * TO O KM BY INSISTING THAT THE VALUE AT O KM FOR A PARTICULAR * PROFILE IS SET EQUAL TO THE CORRESPONDING VALUE AT 2 KM. THERE * ARE 46 LATITUDE DATA POINTS FOR EACH ATMOSPHERIC CONDITION. * LATITUDE COVERAGE IS FROM -90 to +90 AT 4 DEGREE INCREMENTS. * ************************************************************************ ************************************************************************ * * (c) Copyright 1987-2001 by GATS, Inc. * 11864 Canon Blvd. Suite 101 * Newport News, VA 23606 * Phone: (757) 873-5920 * * All Rights Reserved. No part of this software or publication may be * reproduced, stored in a retrieval system, or transmitted, in any form * or by any means, electronic, mechanical, photocopying, recording, or * otherwise without the prior written permission of GATS, Inc. * ************************************************************************ IMPLICIT NONE REAL*4 R,DZ,RI,WGHT,WGHT2,GREF PARAMETER (R = 287.0, GREF=9.8) INTEGER*4 MAXLEV,CO2NUM PARAMETER (MAXLEV = 1000) LOGICAL IVAR,IQUIT INTEGER*4 NTIMES,LEVELS,MR,MGM,IDG(MGM) INTEGER*4 J,IFIX,IPTIE INTEGER*4 NLATITUDES,NLEVELS,I,II,III,IG,II2 REAL*4 G,RE,REARTH,XLAT,CO2I,CO2MAX,ZMBEG,ZMEND,FACT,XLATINC,ZG REAL*4 Q(MR,MGM),ZR(MR),PR(MR),TR(MR),CO2(MR),CO2_FIX,CO2_FACT REAL*4 Z_WACCM(NLEVELS),P_WACCM(NLEVELS) REAL*4 ZGEO_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 T_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 H2O_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 CO2_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 O3_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 N2O_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 CH4_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 O_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 O2_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 NO_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 NO2_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 OH_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 N2_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 HO2_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 CFCL3_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 CF2CL2_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 H_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 O1D_WACCM(NLEVELS,NLATITUDES,NTIMES) REAL*4 EVAR(MAXLEV) REAL*4 CO2PERB(MAXLEV) REAL*4 CO2TEMP(MAXLEV) IF( MR .LE. NLEVELS )THEN WRITE(*,*) 'MR must be at least ',NLEVELS+1 STOP 'ERROR in ATMSTATE' ENDIF IF( MAXLEV .LT. MR )THEN WRITE(*,*) 'MAXLEV must be at least ',MR STOP 'ERROR in ATMSTATE' ENDIF IF( MGM .LT. 16 )THEN WRITE(*,*) 'MGM must be at least 16' STOP 'ERROR in ATMSTATE' ENDIF * * ... Determine index for tangent-point latitude * I=1 IF(XLAT.LT.-90.0.OR.XLAT.GT.90.0) THEN WRITE(*,*) 'ERROR: XLAT not in range' ENDIF XLATINC=180.0/(NLATITUDES-1) RI=1.0+(XLAT+90.0)/XLATINC II=INT(RI) IF(II.LT.1) THEN II=1 WGHT=1.0 ELSEIF(II.GT.NLATITUDES-1) THEN II=NLATITUDES-1 WGHT=0.0 ELSE WGHT=1.0-(RI-REAL(II)) ENDIF IF(WGHT.LT.0.0.OR.WGHT.GT.1.0) THEN WRITE(*,*) 'ERROR: XLAT not in range' ENDIF II2=II+1 WGHT2=1.0-WGHT CTM Set the q values. Not all are available now. DO III = 1, NLEVELS PR(III) = P_WACCM(III) CTM Set ZR at this point to Geopotential altitude ZR(III) = WGHT*ZGEO_WACCM(III,II,I)+WGHT2*ZGEO_WACCM(III,II2,I) TR(III) = WGHT*T_WACCM(III,II,I)+WGHT2*T_WACCM(III,II2,I) Q(III,1) = WGHT*O2_WACCM(III,II,I)+WGHT2*O2_WACCM(III,II2,I) Q(III,2) = WGHT*CO2_WACCM(III,II,I)+WGHT2*CO2_WACCM(III,II2,I) Q(III,3) = WGHT*O3_WACCM(III,II,I)+WGHT2*O3_WACCM(III,II2,I) Q(III,4) = WGHT*H2O_WACCM(III,II,I)+WGHT2*H2O_WACCM(III,II2,I) Q(III,5) = WGHT*N2O_WACCM(III,II,I)+WGHT2*N2O_WACCM(III,II2,I) Q(III,6) = WGHT*CH4_WACCM(III,II,I)+WGHT2*CH4_WACCM(III,II2,I) Q(III,7) = WGHT*NO_WACCM(III,II,I)+WGHT2*NO_WACCM(III,II2,I) Q(III,8) = WGHT*NO2_WACCM(III,II,I)+WGHT2*NO2_WACCM(III,II2,I) Q(III,9) = WGHT*OH_WACCM(III,II,I)+WGHT2*OH_WACCM(III,II2,I) Q(III,10) = WGHT*N2_WACCM(III,II,I)+WGHT2*N2_WACCM(III,II2,I) Q(III,11) = WGHT*HO2_WACCM(III,II,I)+WGHT2*HO2_WACCM(III,II2,I) Q(III,12) = WGHT*O_WACCM(III,II,I)+WGHT2*O_WACCM(III,II2,I) Q(III,13) = WGHT*H_WACCM(III,II,I)+WGHT2*H_WACCM(III,II2,I) Q(III,14) = WGHT*CFCL3_WACCM(III,II,I)+ & WGHT2*CFCL3_WACCM(III,II2,I) Q(III,15 )= WGHT*CF2CL2_WACCM(III,II,I)+ & WGHT2*CF2CL2_WACCM(III,II2,I) Q(III,16) = WGHT*O1D_WACCM(III,II,I)+WGHT2*O1D_WACCM(III,II2,I) ENDDO CALL RADGRAV(0.0,XLAT,G,RE) CO2MAX=0.0 CO2NUM=0 DO III=1,NLEVELS CO2I=Q(III,2) CTM calculate ZR from Geopotential Altitude, 0.001 factor converts to km ZG=0.001*ZR(III)*GREF/G C write(99,*) III,RE,ZG ZR(III)=RE*ZG/(RE-ZG) C write(99,*) III,ZR(III) IF(PR(III).GT.1.0.AND.PR(III).LT.10.0) THEN CO2MAX=CO2MAX+CO2I CO2NUM=CO2NUM+1 ENDIF ENDDO CO2MAX=CO2MAX/CO2NUM WRITE(*,*) 'CO2: TREND, WACCM_REFC2',CO2_FIX,CO2MAX CTM USE WACCM as is C FIX CO2 PROFILE CTM First check whether CO2 VMR exceeds CO2_FIX. If so then scale upper CTM altitudes using last value .LE. CO2_FIX C CO2MAX=0 C DO III=1,NLEVELS C CO2I=Q(III,2) C IF(PR(III).GT.0.05.AND.PR(III).LT.0.5) THEN C CO2MAX=MAX(CO2MAX,CO2I) C IF(CO2I.GT.CO2_FIX) THEN C CO2_FACT=CO2_FIX/CO2I C IFIX=III C GOTO 11 C ENDIF C ENDIF C ENDDO CTM IF the above condition isn't met then scale upper altitudes using CTM maximum value in array as long as it's close enough C CO2_FACT=CO2_FIX/CO2MAX C IF(CO2_FACT.LT.1.02) THEN C DO III=1,NLEVELS C CO2I=Q(III,2) C IF(PR(III).GT.0.05.AND.PR(III).LT.0.5.AND.CO2I.EQ.CO2MAX) THEN C IFIX=III C GOTO 11 C ENDIF C ENDDO C ENDIF C WRITE(*,*) 'ERROR FIXING CO2 PROFILE' C STOP 'ERROR IN ATMSTATE' C 11 DO III=1,NLEVELS C IF(III.LT.IFIX) THEN C Q(III,2)=CO2_FACT*Q(III,2) C ELSE C Q(III,2)=CO2_FIX C ENDIF C CO2(III)=Q(III,2) C ENDDO CTM Perturb CO2 profile above 70km. Start with 0% change at 70km and then CTM slowly increase to 15% at 80km and keep it at 15% above 80km. CTM Do this by first creating CO2PERB=CO2*1.15 and splicing with CO2 CTM using ZMBEG=80km and ZMEND=70km. Write out warning to standard out CTM so it is obvious this is being done. CTM Set up parameters for splicing perturbed CO2 to WACCM for test CTM TURN THIS OFF!!! C IVAR = .FALSE. C ZMBEG = 80.0 C ZMEND = 70.0 C FACT=1.0 C WRITE(*,*) "!!!!WARNING: Perturbing CO2 for test run :WARNING!!!!" C DO III=1,NLEVELS C CO2PERB(III)=CO2(III)*1.15 C ENDDO C CALL FUSE(MR,NLEVELS,ZMBEG,ZMEND,ZR,CO2,CO2PERB,FACT,IVAR, C & EVAR,EVAR,CO2TEMP,IQUIT) C DO III=1,NLEVELS C WRITE(*,*) ZR(III),CO2PERB(III),CO2(III),CO2TEMP(III) C CO2(III)=CO2TEMP(III) C ENDDO CTM END CO2 Perturb * ZR(MR) = 0.0 DZ = 1000.0*ZR(NLEVELS) TR(MR) = TR(NLEVELS) PR(MR) = PR(NLEVELS)*EXP(G*DZ/R/TR(NLEVELS)) CO2(MR)= CO2(NLEVELS) DO IG = 1, MGM Q(MR,IG) = Q(NLEVELS,IG) ENDDO CTM Do not call Hypsometric, Pressures are not hydrostatic at this point. CTM HYPSOMETRIC recalculates altitudes for hydrostatic balance CTM Hysometric integrates from top need to integrate from bottom C CALL HYPSOMETRIC(1,MR,MR,RE,XLAT,MR,MGM,IDG,ZR,PR,TR,Q) C IPTIE=MR C DO III=1,MR C IF(PR(III).LT.10.0) IPTIE=III C ENDDO C CALL HYDROZ(1,MR,IPTIE,RE,XLAT,ZR,PR,TR) * * DO III = 1, MR * WRITE(45,101) III, II, IREC, DAY, ZR(III), PR(III), TR(III) * ENDDO * *101 FORMAT(I3,2X,I2,2X,I3,2X,I3,2X,1PE12.5,2X,1PE12.5,2X,1PE12.5) * DO I = 1, MR * WRITE(46,102) I,ZR(I),Q(I,1),Q(I,2),Q(I,3),Q(I,4) * ENDDO *102 FORMAT(I2,5(1X,1PE12.5)) * RETURN END