!-- 130204,add Time /= 99.99 check in chkinp() for daily 1deg processing. !-- 120425,: ! When compile, use pgf95. Dont use ifort120, it'll generate a ! buggy code that product wired invalide Dtauex values. !-- 160410, : added capability to allow uer provided aerosol model. !-- 200628, : Modified the use of "SRB_UserAer_Form1" or "SRB_UserAer_Form2" type; ! Wrapped all subroutine in a Fortran Module so interface checking is enforced; ! Removed implcit type variables. ! Note: "McClatchey_Atmospheres" and "Standard_Aerosol" modules have to be compiled ! before the main module "Module_SRB", so they are placed infront of the "Module_SRB". !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- MODULE McClatchey_Atmospheres !------------------------------------------------------------------------------- C .. Implicit None Statement .. IMPLICIT NONE C .. C .. Parameters .. INTEGER Mccl PARAMETER ( Mccl=31 ) C .. C .. Local Scalars .. INTEGER, PRIVATE :: L C .. C .. Local Arrays .. CHARACTER Atm( 6 ) * 18 REAL H2o( Mccl, 6 ), O3( Mccl, 6 ), Ps( Mccl, 6 ), Ts( Mccl, 6 ), & Zs( Mccl ) C .. C .. Data statements .. C *** MCCLATCHEY ALTITUDES DATA Zs / 50., 45., 40., 35., 30., 25., 24., 23., 22., 21., 20., & 19., 18., 17., 16., 15., 14., 13., 12., 11., 10., 9., 8., & 7., 6., 5., 4., 3., 2., 1., 0. / C *** T R O P I C A L DATA ( Ps(L,1), L=1, Mccl ) / 0.854, 1.59, 3.05, 6.00, 12.2, & 25.7, 30.0, 35.0, 40.9, 48.0, 56.5, 66.6, 78.9, 93.7, 111.0, & 132.0, 156.0, 182.0, 213.0, 247.0, 286.0, 329.0, 378.0, & 432.0, 492.0, 559.0, 633.0, 715.0, 805.0, 904.0, 1013. / DATA ( Ts(L,1), L=1, Mccl ) / 270., 265., 254., 243., 232., 221., & 219., 217., 215., 211., 207., 203., 199., 195., 197., 204., & 210., 217., 224., 230., 237., 244., 250., 257., 264., 270., & 277., 284., 288., 294., 300. / DATA ( H2o(L,1), L=1, Mccl ) / 6.3E-6, 1.9E-5, 4.3E-5, 1.1E-4, & 3.6E-4, 6.7E-4, 6.0E-4, 5.4E-4, 5.1E-4, 5.1E-4, 4.5E-4, & 4.9E-4, 5.0E-4, 5.6E-4, 6.4E-4, 7.6E-4, 1.0E-3, 1.8E-3, & 6.0E-3, 1.7E-2, 5.0E-2, 0.12, 0.25, 0.47, 0.85, 1.5, 2.2, & 4.7, 9.3, 13., 19. / DATA ( O3(L,1), L=1, Mccl ) / 4.3E-6, 1.3E-5, 4.1E-5, 9.2E-5, & 2.4E-4, 3.4E-4, 3.4E-4, 3.2E-4, 2.8E-4, 2.4E-4, 1.9E-4, & 1.4E-4, 9.0E-5, 6.9E-5, 4.7E-5, 4.7E-5, 2 * 4.5E-5, 4.3E-5, & 4.1E-5, 3 * 3.9E-5, 4.1E-5, 4.3E-5, 4.5E-5, 4.7E-5, 5.1E-5, & 5.4E-5, 2 * 5.6E-5 / C *** M I D L A T I T U D E S U M M E R DATA ( Ps(L,2), L=1, Mccl ) / 0.951, 1.76, 3.33, 6.52, 13.2, & 27.7, 32.2, 37.6, 43.7, 51., 59.5, 69.5, 81.2, 95., 111., & 130., 153., 179., 209., 243., 281., 324., 372., 426., 487., & 554., 628., 710., 802., 902., 1013. / DATA ( Ts(L,2), L=1, Mccl ) / 276., 270., 258., 245., 234., 224., & 223., 222., 220., 219., 218., 217., 6 * 216., 222., 229., & 235., 242., 248., 255., 261., 267., 273., 279., 285., 290., & 294. / DATA ( H2o(L,2), L=1, Mccl ) / 6.3E-6, 1.9E-5, 4.3E-5, 1.1E-4, & 3.6E-4, 6.7E-4, 6.0E-4, 5.4E-4, 5.1E-4, 5.1E-4, 4.5E-4, & 4.9E-4, 5.0E-4, 5.6E-4, 6.4E-4, 7.6E-4, 1.0E-3, 1.8E-3, & 6.0E-3, 2.2E-2, 6.4E-2, 0.12, 0.21, 0.37, 0.61, 1., 1.9, & 3.3, 5.9, 9.3, 14. / DATA ( O3(L,2), L=1, Mccl ) / 4.3E-6, 1.3E-5, 4.1E-5, 9.2E-5, & 2.0E-4, 3.0E-4, 3.2E-4, 3.4E-4, 3.6E-4, 3.6E-4, 3.4E-4, & 3.2E-4, 2.8E-4, 2.4E-4, 2.1E-4, 1.9E-4, 1.8E-4, 1.5E-4, & 1.2E-4, 1.1E-4, 9.0E-5, 8.6E-5, 7.9E-5, 7.5E-5, 6.9E-5, & 6.6E-5, 6.4E-5, 6.2E-5, 3 * 6.0E-5 / C *** M I D L A T I T U D E W I N T E R DATA ( Ps(L,3), L=1, Mccl ) / 0.682, 1.29, 2.53, 5.18, 11.1, & 24.3, 28.6, 33.4, 39.1, 45.8, 53.7, 62.8, 73.5, 86.1, 100.7, & 117.8, 137.8, 161., 188.2, 219.9, 256.8, 299.2, 347.3, & 401.6, 462.7, 531.3, 608.1, 693.8, 789.7, 897.3, 1018. / DATA ( Ts(L,3), L=1, Mccl ) / 265.7, 258.5, 243.2, 227.8, 217.4, & 7 * 215.2, 215.7, 216.2, 216.7, 217.2, 217.7, 218.2, 218.7, & 219.2, 219.7, 225.7, 231.7, 237.7, 243.7, 249.7, 255.7, & 261.7, 265.2, 268.7, 272.2 / DATA ( H2o(L,3), L=1, Mccl ) / 6.3E-6, 1.9E-5, 4.3E-5, 1.1E-4, & 3.6E-4, 6.7E-4, 6.0E-4, 5.4E-4, 5.1E-4, 5.1E-4, 4.5E-4, & 4.9E-4, 5.0E-4, 5.6E-4, 6.4E-4, 7.6E-4, 1.0E-3, 1.8E-3, & 6.0E-3, 6.9E-3, 7.5E-3, 1.6E-2, 3.5E-2, 8.5E-2, 0.21, 0.38, & 0.66, 1.2, 1.8, 2.5, 3.5 / DATA ( O3(L,3), L=1, Mccl ) / 4.3E-6, 1.3E-5, 4.1E-5, 9.2E-5, & 1.9E-4, 3.4E-4, 3.6E-4, 3.9E-4, 2 * 4.3E-4, 4.5E-4, 4.3E-4, & 4.1E-4, 3.9E-4, 3.6E-4, 3.4E-4, 3.2E-4, 3.0E-4, 2.6E-4, & 2.1E-4, 1.6E-4, 1.2E-4, 9.0E-5, 7.7E-5, 6.4E-5, 5.8E-5, & 3 * 4.9E-5, 5.4E-5, 6.0E-5 / C *** S U B A R C T I C S U M M E R DATA ( Ps(L,4), L=1, Mccl ) / 0.987, 1.81, 3.4, 6.61, 13.4, 27.8, & 32.27, 37.5, 43.6, 50.7, 58.9, 68.6, 79.8, 92.8, 108., 125., & 146., 170., 197.7, 230., 267.7, 310.7, 359., 413., 473., & 541., 616., 700., 792.9, 896., 1010. / DATA ( Ts(L,4), L=1, Mccl ) / 277., 274., 262., 247., 235., 228., & 226., 14 * 225., 232., 239., 246., 253., 260., 266., 271., & 276., 282., 287. / DATA ( H2o(L,4), L=1, Mccl ) / 6.3E-6, 1.9E-5, 4.3E-5, 1.1E-4, & 3.6E-4, 6.7E-4, 6.0E-4, 5.4E-4, 5.1E-4, 5.1E-4, 4.5E-4, & 4.9E-4, 5.0E-4, 5.6E-4, 6.4E-4, 7.6E-4, 1.0E-3, 1.8E-3, & 6.0E-3, 9.4E-3, 1.5E-2, 4.2E-2, 0.13, 0.29, 0.54, 1.0, 1.7, & 2.7, 4.2, 6.0, 9.1 / DATA ( O3(L,4), L=1, Mccl ) / 4.3E-6, 1.3E-5, 4.1E-5, 9.2E-5, & 1.4E-4, 2.6E-4, 2.8E-4, 3.0E-4, 3.2E-4, 3.6E-4, 3.9E-4, & 2 * 4.1E-4, 3.9E-4, 3.4E-4, 3.2E-4, 2.8E-4, 2.6E-4, 2.1E-4, & 1.8E-4, 1.3E-4, 1.1E-4, 7.9E-5, 7.5E-5, 7.1E-5, 6.4E-5, & 6.0E-5, 5.8E-5, 5.6E-5, 5.4E-5, 4.9E-5 / C *** S U B A R C T I C W I N T E R DATA ( Ps(L,5), L=1, Mccl ) / 0.5719, 1.113, 2.243, 4.701, 10.2, & 22.56, 26.49, 31.09, 36.47, 42.77, 50.14, 58.75, 68.82, & 80.58, 94.31, 110.3, 129.1, 151., 176.6, 206.7, 241.8, & 282.9, 330.8, 385.3, 446.7, 515.8, 593.2, 679.8, 777.5, & 887.8, 1013. / DATA ( Ts(L,5), L=1, Mccl ) / 259.3, 247., 234.7, 222.2, 216., & 211.2, 211.8, 212.4, 213., 213.6, 214.1, 214.8, 215.4, 216., & 216.6, 7 * 217.2, 220.6, 227.3, 234.1, 240.9, 247.7, 252.7, & 255.9, 259.1, 257.1 / DATA ( H2o(L,5), L=1, Mccl ) / 6.3E-6, 1.9E-5, 4.3E-5, 1.1E-4, & 3.6E-4, 6.7E-4, 6.0E-4, 5.4E-4, 5.1E-4, 5.1E-4, 4.5E-4, & 4.9E-4, 5.0E-4, 5.6E-4, 6.4E-4, 7.6E-4, 1.0E-3, 1.8E-3, & 2.6E-3, 3.8E-3, 5.5E-3, 8.4E-3, 1.1E-2, 5.4E-2, 9.8E-2, & 0.20, 0.41, 0.68, 0.94, 1.2, 1.2 / DATA ( O3(L,5), L=1, Mccl ) / 4.3E-6, 1.3E-5, 4.1E-5, 9.2E-5, & 1.5E-4, 3.2E-4, 3.6E-4, 4.3E-4, 4.7E-4, 5.1E-4, 5.6E-4, & 6.0E-4, 3 * 6.2E-4, 5.6E-4, 4.9E-4, 4.7E-4, 4.3E-4, 3.2E-4, & 2.4E-4, 1.6E-4, 9.0E-5, 7.1E-5, 4.9E-5, 4.7E-5, 4.5E-5, & 4.3E-5, 3 * 4.1E-5 / C *** U. S. S T A N D A R D DATA ( Ps(L,6), L=1, Mccl ) / 0.798, 1.491, 2.871, 5.746, 11.97, & 25.49, 29.72, 34.67, 40.47, 47.29, 55.29, 64.67, 75.65, & 88.5, 103.5, 121.1, 141.7, 165.8, 194., 227., 265., 308., & 356.5, 411.1, 472.2, 540.5, 616.6, 701.2, 795., 898.6, & 1013. / DATA ( Ts(L,6), L=1, Mccl ) / 270.6, 264.2, 253.4, 236.5, 226.5, & 221.6, 220.6, 219.6, 218.6, 217.6, 9 * 216.6, 216.8, 223.2, & 229.7, 236.2, 242.7, 249.2, 255.7, 262.2, 268.7, 275.1, & 281.6, 288.1 / DATA ( H2o(L,6), L=1, Mccl ) / 1.2E-5, 3.2E-5, 6.7E-5, 1.6E-4, & 3.8E-4, 6.6E-4, 6.1E-4, 5.7E-4, 5.2E-4, 4.8E-4, 3 * 4.4E-4, & 5.2E-4, 6.1E-4, 7.2E-4, 8.4E-4, 1.8E-3, 3.7E-3, 8.2E-3, & 1.8E-2, 4.6E-2, 0.12, 0.21, 0.38, 0.64, 1.1, 1.8, 2.9, 4.2, & 5.9 / DATA ( O3(L,6), L=1, Mccl ) / 4.0E-6, 1.7E-5, 4.9E-5, 1.1E-4, & 2.0E-4, 3.4E-4, 3.6E-4, 3.8E-4, 3.9E-4, 2 * 3.8E-4, 3.5E-4, & 3.2E-4, 2.8E-4, 2.4E-4, 2.1E-4, 1.9E-4, 1.7E-4, 1.6E-4, & 1.3E-4, 9.0E-5, 7.1E-5, 5.2E-5, 4.9E-5, 4.5E-5, 2 * 4.6E-5, & 5.0E-5, 3 * 5.4E-5 / DATA Atm / 'TROPICAL ', 'MIDLATITUDE SUMMER', & 'MIDLATITUDE WINTER', 'SUBARCTIC SUMMER ', & 'SUBARCTIC WINTER ', 'U.S. STANDARD ' / END MODULE McClatchey_Atmospheres !----------------------------------------------------------------------- !----------------------------------------------------------------------- MODULE Standard_Aerosol !----------------------------------------------------------------------- C .. Implicit None Statement .. IMPLICIT NONE PRIVATE PUBLIC :: Mlsra, Nzsra, Za, & Bsra, Gsra, Osra, Sclht C .. C .. Parameters .. C Mlsra : max number of layers in standard radiation atmospheres C Mzsra : max number of levels in standard radiation atmospheres C Mnwl : max number of spectral intervals in aerosol models C Msra : max number of aerosol models INTEGER Msra, Mzsra, Mlsra, Mnwl PARAMETER ( Msra=5, Mzsra=4, Mlsra=Mzsra-1, Mnwl=7 ) C .. C .. Local Scalars .. INTEGER, PRIVATE :: L C .. C .. Local Arrays .. CHARACTER Aerosl( Msra ) * 70 INTEGER Nzsra( Msra ) REAL Bsra( Mnwl, Mlsra, Msra ), Gsra( Mnwl, Mlsra, Msra ), & Osra( Mnwl, Mlsra, Msra ), Za( Mzsra, Msra ), & Sclht( Mlsra, Msra ) C .. C .. Data statements .. DATA Nzsra / 4, 4, 4, 4, 4 / DATA Aerosl( 1 ) / '*** AEROSOL PROFILE CONT-I' / DATA ( Za(L,1), L=1, 4 ) / 35., 12., 2., 0. / DATA ( Sclht(L,1), L=1, 3 ) / 99., 8., 8. / DATA ( Bsra(1,L,1), L=1, 3 ) / 2.99E-4, 5.16E-3, 1.07E-1 / DATA ( Bsra(2,L,1), L=1, 3 ) / 2.54E-4, 3.90E-3, 8.19E-2 / DATA ( Bsra(3,L,1), L=1, 3 ) / 2.14E-4, 2.98E-3, 6.38E-2 / DATA ( Bsra(4,L,1), L=1, 3 ) / 1.80E-4, 2.33E-3, 5.08E-2 / DATA ( Bsra(5,L,1), L=1, 3 ) / 1.19E-4, 1.40E-3, 3.22E-2 / DATA ( Bsra(6,L,1), L=1, 3 ) / 4.31E-5, 4.80E-4, 1.36E-2 / DATA ( Bsra(7,L,1), L=1, 3 ) / 1.64E-5, 1.82E-4, 6.84E-3 / DATA ( Osra(1,L,1), L=1, 3 ) / 1.000, 0.901, 0.917 / DATA ( Osra(2,L,1), L=1, 3 ) / 1.000, 0.901, 0.918 / DATA ( Osra(3,L,1), L=1, 3 ) / 1.000, 0.894, 0.912 / DATA ( Osra(4,L,1), L=1, 3 ) / 1.000, 0.885, 0.905 / DATA ( Osra(5,L,1), L=1, 3 ) / 0.999, 0.851, 0.875 / DATA ( Osra(6,L,1), L=1, 3 ) / 0.997, 0.777, 0.805 / DATA ( Osra(7,L,1), L=1, 3 ) / 0.408, 0.504, 0.525 / DATA ( Gsra(1,L,1), L=1, 3 ) / 0.717, 0.699, 0.715 / DATA ( Gsra(2,L,1), L=1, 3 ) / 0.723, 0.686, 0.703 / DATA ( Gsra(3,L,1), L=1, 3 ) / 0.717, 0.674, 0.690 / DATA ( Gsra(4,L,1), L=1, 3 ) / 0.707, 0.661, 0.678 / DATA ( Gsra(5,L,1), L=1, 3 ) / 0.683, 0.641, 0.655 / DATA ( Gsra(6,L,1), L=1, 3 ) / 0.620, 0.633, 0.636 / DATA ( Gsra(7,L,1), L=1, 3 ) / 0.485, 0.758, 0.729 / DATA Aerosl( 2 ) / '*** AEROSOL PROFILE DESERT' / DATA ( Za(L,2), L=1, 4 ) / 35., 12., 6., 0. / DATA ( Sclht(L,2), L=1, 3 ) / 99., 8., 2. / DATA ( Bsra(1,L,2), L=1, 3 ) / 2.99E-4, 5.16E-3, 1.53E-1 / DATA ( Bsra(2,L,2), L=1, 3 ) / 2.54E-4, 3.90E-3, 1.46E-1 / DATA ( Bsra(3,L,2), L=1, 3 ) / 2.14E-4, 2.98E-3, 1.41E-1 / DATA ( Bsra(4,L,2), L=1, 3 ) / 1.80E-4, 2.33E-3, 1.38E-1 / DATA ( Bsra(5,L,2), L=1, 3 ) / 1.19E-4, 1.40E-3, 1.33E-1 / DATA ( Bsra(6,L,2), L=1, 3 ) / 4.31E-5, 4.80E-4, 1.22E-1 / DATA ( Bsra(7,L,2), L=1, 3 ) / 1.64E-5, 1.82E-4, 9.28E-2 / DATA ( Osra(1,L,2), L=1, 3 ) / 1.000, 0.901, 0.761 / DATA ( Osra(2,L,2), L=1, 3 ) / 1.000, 0.901, 0.831 / DATA ( Osra(3,L,2), L=1, 3 ) / 1.000, 0.894, 0.881 / DATA ( Osra(4,L,2), L=1, 3 ) / 1.000, 0.885, 0.908 / DATA ( Osra(5,L,2), L=1, 3 ) / 0.999, 0.851, 0.926 / DATA ( Osra(6,L,2), L=1, 3 ) / 0.997, 0.777, 0.928 / DATA ( Osra(7,L,2), L=1, 3 ) / 0.408, 0.504, 0.857 / DATA ( Gsra(1,L,2), L=1, 3 ) / 0.717, 0.699, 0.785 / DATA ( Gsra(2,L,2), L=1, 3 ) / 0.723, 0.686, 0.753 / DATA ( Gsra(3,L,2), L=1, 3 ) / 0.717, 0.674, 0.729 / DATA ( Gsra(4,L,2), L=1, 3 ) / 0.707, 0.661, 0.713 / DATA ( Gsra(5,L,2), L=1, 3 ) / 0.683, 0.641, 0.698 / DATA ( Gsra(6,L,2), L=1, 3 ) / 0.620, 0.633, 0.691 / DATA ( Gsra(7,L,2), L=1, 3 ) / 0.485, 0.758, 0.701 / DATA Aerosl( 3 ) / '*** AEROSOL PROFILE MAR-I' / DATA ( Za(L,3), L=1, 4 ) / 35., 12., 2., 0. / DATA ( Sclht(L,3), L=1, 3 ) / 99., 8., 1. / DATA ( Bsra(1,L,3), L=1, 3 ) / 2.99E-4, 5.16E-3, 7.94E-2 / DATA ( Bsra(2,L,3), L=1, 3 ) / 2.54E-4, 3.90E-3, 7.68E-2 / DATA ( Bsra(3,L,3), L=1, 3 ) / 2.14E-4, 2.98E-3, 7.52E-2 / DATA ( Bsra(4,L,3), L=1, 3 ) / 1.80E-4, 2.33E-3, 7.40E-2 / DATA ( Bsra(5,L,3), L=1, 3 ) / 1.19E-4, 1.40E-3, 7.10E-2 / DATA ( Bsra(6,L,3), L=1, 3 ) / 4.31E-5, 4.80E-4, 5.82E-2 / DATA ( Bsra(7,L,3), L=1, 3 ) / 1.64E-5, 1.82E-4, 3.84E-2 / DATA ( Osra(1,L,3), L=1, 3 ) / 1.000, 0.901, 0.995 / DATA ( Osra(2,L,3), L=1, 3 ) / 1.000, 0.901, 0.997 / DATA ( Osra(3,L,3), L=1, 3 ) / 1.000, 0.894, 0.997 / DATA ( Osra(4,L,3), L=1, 3 ) / 1.000, 0.885, 0.997 / DATA ( Osra(5,L,3), L=1, 3 ) / 0.999, 0.851, 0.996 / DATA ( Osra(6,L,3), L=1, 3 ) / 0.997, 0.777, 0.993 / DATA ( Osra(7,L,3), L=1, 3 ) / 0.408, 0.504, 0.791 / DATA ( Gsra(1,L,3), L=1, 3 ) / 0.717, 0.699, 0.769 / DATA ( Gsra(2,L,3), L=1, 3 ) / 0.723, 0.686, 0.766 / DATA ( Gsra(3,L,3), L=1, 3 ) / 0.717, 0.674, 0.765 / DATA ( Gsra(4,L,3), L=1, 3 ) / 0.707, 0.661, 0.767 / DATA ( Gsra(5,L,3), L=1, 3 ) / 0.683, 0.641, 0.774 / DATA ( Gsra(6,L,3), L=1, 3 ) / 0.620, 0.633, 0.787 / DATA ( Gsra(7,L,3), L=1, 3 ) / 0.485, 0.758, 0.775 / DATA Aerosl( 4 ) / '*** AEROSOL PROFILE ARCTIC' / DATA ( Za(L,4), L=1, 4 ) / 35., 12., 2., 0. / DATA ( Sclht(L,4), L=1, 3 ) / 99., 8., 99. / DATA ( Bsra(1,L,4), L=1, 3 ) / 2.99E-4, 5.16E-3, 2.86E-2 / DATA ( Bsra(2,L,4), L=1, 3 ) / 2.54E-4, 3.90E-3, 2.32E-2 / DATA ( Bsra(3,L,4), L=1, 3 ) / 2.14E-4, 2.98E-3, 1.94E-2 / DATA ( Bsra(4,L,4), L=1, 3 ) / 1.80E-4, 2.33E-3, 1.66E-2 / DATA ( Bsra(5,L,4), L=1, 3 ) / 1.19E-4, 1.40E-3, 1.27E-2 / DATA ( Bsra(6,L,4), L=1, 3 ) / 4.31E-5, 4.80E-4, 7.80E-3 / DATA ( Bsra(7,L,4), L=1, 3 ) / 1.64E-5, 1.82E-4, 4.66E-3 / DATA ( Osra(1,L,4), L=1, 3 ) / 1.000, 0.901, 0.854 / DATA ( Osra(2,L,4), L=1, 3 ) / 1.000, 0.901, 0.863 / DATA ( Osra(3,L,4), L=1, 3 ) / 1.000, 0.894, 0.868 / DATA ( Osra(4,L,4), L=1, 3 ) / 1.000, 0.885, 0.872 / DATA ( Osra(5,L,4), L=1, 3 ) / 0.999, 0.851, 0.874 / DATA ( Osra(6,L,4), L=1, 3 ) / 0.997, 0.777, 0.876 / DATA ( Osra(7,L,4), L=1, 3 ) / 0.408, 0.504, 0.667 / DATA ( Gsra(1,L,4), L=1, 3 ) / 0.717, 0.699, 0.716 / DATA ( Gsra(2,L,4), L=1, 3 ) / 0.723, 0.686, 0.711 / DATA ( Gsra(3,L,4), L=1, 3 ) / 0.717, 0.674, 0.709 / DATA ( Gsra(4,L,4), L=1, 3 ) / 0.707, 0.661, 0.709 / DATA ( Gsra(5,L,4), L=1, 3 ) / 0.683, 0.641, 0.717 / DATA ( Gsra(6,L,4), L=1, 3 ) / 0.620, 0.633, 0.743 / DATA ( Gsra(7,L,4), L=1, 3 ) / 0.485, 0.758, 0.752 / DATA Aerosl( 5 ) / '*** AEROSOL PROFILE ANTARCTIC' / DATA ( Za(L,5), L=1, 4 ) / 35., 12., 10., 0. / DATA ( Sclht(L,5), L=1, 3 ) / 99., 8., 8. / DATA ( Bsra(1,L,5), L=1, 3 ) / 2.99E-4, 5.16E-3, 1.14E-2 / DATA ( Bsra(2,L,5), L=1, 3 ) / 2.54E-4, 3.90E-3, 1.05E-2 / DATA ( Bsra(3,L,5), L=1, 3 ) / 2.14E-4, 2.98E-3, 9.45E-3 / DATA ( Bsra(4,L,5), L=1, 3 ) / 1.80E-4, 2.33E-3, 8.39E-3 / DATA ( Bsra(5,L,5), L=1, 3 ) / 1.19E-4, 1.40E-3, 6.17E-3 / DATA ( Bsra(6,L,5), L=1, 3 ) / 4.31E-5, 4.80E-4, 2.81E-3 / DATA ( Bsra(7,L,5), L=1, 3 ) / 1.64E-5, 1.82E-4, 1.32E-3 / DATA ( Osra(1,L,5), L=1, 3 ) / 1.000, 0.901, 0.999 / DATA ( Osra(2,L,5), L=1, 3 ) / 1.000, 0.901, 0.999 / DATA ( Osra(3,L,5), L=1, 3 ) / 1.000, 0.894, 0.999 / DATA ( Osra(4,L,5), L=1, 3 ) / 1.000, 0.885, 0.999 / DATA ( Osra(5,L,5), L=1, 3 ) / 0.999, 0.851, 0.999 / DATA ( Osra(6,L,5), L=1, 3 ) / 0.997, 0.777, 0.997 / DATA ( Osra(7,L,5), L=1, 3 ) / 0.408, 0.504, 0.544 / DATA ( Gsra(1,L,5), L=1, 3 ) / 0.717, 0.699, 0.774 / DATA ( Gsra(2,L,5), L=1, 3 ) / 0.723, 0.686, 0.779 / DATA ( Gsra(3,L,5), L=1, 3 ) / 0.717, 0.674, 0.778 / DATA ( Gsra(4,L,5), L=1, 3 ) / 0.707, 0.661, 0.774 / DATA ( Gsra(5,L,5), L=1, 3 ) / 0.683, 0.641, 0.759 / DATA ( Gsra(6,L,5), L=1, 3 ) / 0.620, 0.633, 0.715 / DATA ( Gsra(7,L,5), L=1, 3 ) / 0.485, 0.758, 0.608 / END MODULE Standard_Aerosol !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- MODULE Module_SRB !------------------------------------------------------------------------------- IMPLICIT NONE PRIVATE PUBLIC :: & Mnwl, Mlsra, Mzsra, & SRB_UserAer_Form1, & SRB_UserAer_Form2, & Modflx !--- User aerosol input form 1 ! integer ,PARAMETER :: Mnwl=7 !max number of spectral intervals in aerosol models integer ,PARAMETER :: Mlsra=3 !max number of aer layers. integer ,PARAMETER :: Mzsra=mlsra+1 !max number of aer levels TYPE SRB_UserAer_Form1 real :: Tauaer = 0. !Total column aerosol optical depth character(70) :: AerID = '' !text ID for the aer model integer :: Nzsra = 0 !number of aer profle Z levels real :: Za( Mzsra ) !aer profile level hight (km) real :: Sclht( Mlsra ) !scale height REAL :: Bsra( Mnwl, Mlsra ) !extinction coefficent real :: Osra( Mnwl, Mlsra ) !single scattering albedo real :: Gsra( Mnwl, Mlsra ) !asymmetry factor END TYPE !--- User aerosol input form 2 ! TYPE SRB_UserAer_Form2 real :: Tauaer = 0. !Total column aerosol optical depth real :: absTau = 0. !Total column absorption optical depth character(70) :: AerID = '' !text ID for the aer model END TYPE CONTAINS !------------------------ CONTAINS ------------------------------ !----------------------------------------------------------------------- !----------------------------------------------------------------------- SUBROUTINE Modflx( Prnt, Header, Nwvl, Solcon, Smumin, Smumax, & Misval, Glat, Glon, Month, Time, Jday, & Cwcov,Cicov, Pwater, Ozone, Solamu, Zsrf, & Psrf, aerObj, Tauwcd, Tauicd, Cldwer, Cldier, & Pcldt, Albdir, Albdif, flxall,flxclr ) !!yma 160410 !----------------------------------------------------------------------- C C C Compute the instantaneous shortwave up- and downward fluxes C at specified pressure levels at a single C location from given column-amount of ozone and water vapor, C aerosol and cloud optical properties, spectral surface albedo, and C solar zenith angle. C C Principle of Operation: C C Modflx was designed to provide an efficient way of calculating C solar fluxes from a set of input surface and atomsperic parameters C derived from MODIS instrument on board Terra and Aqua. With small C modification, it can be implemented with parameters from other C satellites. C C INPUT VARIABLES: C C Albdif( Iw ) : Iw=1 to Nwvl; diffuse spectral surface albedo (see C NOTE-1 below) C Albdir( Iw ) : Iw=1 to Nwvl; direct spectral surface albedo (see C NOTE-1 below) C Header : 256-character long description (optional) C Prnt : if TRUE, print input data and fluxes to logical C unit 6. (It is only intended for debugging. C WARNING: for long runs this could produce a C large file, so normally it should be set to C FALSE, and fluxes should be saved in the C calling routine ) C Jday : day of year (1-366) C Nwvl : number of spectral intervals for surface albedo; C in current version, it MUST be 7 (Seven). C Cwcov : water cloud cover fraction (valid range: 0-1) C Cicov : ice cloud cover fraction (valid range: 0-1) C Cldwer : cloud droplet effective radius for water cloud(Micron) C Cldier : cloud particle effective radius for ice cloud(Micron) C Glat : latitude in degrees (valid range: -90 - +90, north C is positive ) C Glon : longitude in degrees (valid range: -180 - +180, C east is positive ) C Misval : special value representing missing atmospheric and C surface (ozone, water vapor, optical depth, and C surface albedo) data C Month : month of year C Ozone : column amount of ozone (atm-cm) C Pwater : column amount of water vapor (g/cm^2) C Pcldt : cloud top pressure C Psrf : surface pressure(mb) C Smumax : maximum solar zenith angle cosine to be used; data C with solar zenith angle cosines larger than C Smumax are skipped C Smumin : minimum solar zenith angle cosine to be used; data C with solar zenith angle cosines smaller than C Smumin are skipped C Solamu : cosine of solar zenith angle at location Glat,Glon C at time Time on day Jday C Solcon : solar constant (W/m^2) for the total spectrum C Tauaer : aerosol optical depth at 0.55 microns C Tauwcd : water cloud optical depth at 0.55 microns C Tauicd : ice cloud optical depth at 0.55 microns C Time : time of atmospheric and surface data (hours) C Zsrf : surface elevation(m) C NOTE-1: surface albedos should be provided for the following C five spectral intervals: (1) 0.2-0.4, (2) 0.4-0.5, C (3) 0.5-0.6, (4) 0.6-0.7, (5) 0.7-1.19, C (6) 1.19-2.38, (7) 2.38-4.00 microns) C C OUTPUT VARIABLES: C C Flxall(Ip,Iv) : Ip=1 to Nfpar, Iv=1 to Nintv, all-sky flux; C Ip=1, top of atmosphere down-flux C Ip=2, top of atmosphere up-flux C Ip=3, surface down-flux C Ip=4, surface up-flux C Ip=5, diffuse surface down-flux C Iv=1, shortwave C Iv=2, visible (PAR) C Iv=3, near infrared C Flxclr(Ip,Iv) : as Flxall but clear-sky flux C C NOTE-2: All fluxes are in W/m^2, C SW = shortwave (0.2-4.0 microns), C PAR = photosynthetically active radiation (0.4-0.7 C microns) C C Data Files: C C In addition to the input variables listed above, the program reads C the following two files from logical units 7 and 8, respectively, in C subroutines Getmod and Setcth: C C vegtyp.dat (vegetation type data of Matthews), C cld_thk.data (Cloud layer thickness model coeffcients). C C These files are expected to be in folder (directory) Auxdata. C C L O C A L V A R I A B L E S: C C C Declin : declination of sun (degrees) C C Flx(Ip,In,Is) : Ip=1 to Nfpar, In=1 Nintv, Is=1 to Nscene; C flux retrived for all scene types (W/m^2) C Fnlst : path and name of output file for listing (and debug) C Fnct : path and name of file of cloud layer thickness model C coefficients C Fnvt : path and name of file containing surface type data C Iaer : aerosol model id (1=Cont; 2=Desert; 3=Mari; C 4=Arctic; 5=Antarctic ) C Iatm : atmospheric model; one of six standard McClatchey C atmospheres: 1=Tropical, 2= Midlatitude summer, C 3=Midlatitude winter, 4=Subarctic summer, C 5=Subarctic winter, 6=US standard C Isrtyp : surface type id (1-water, 2-vegetation, 3-desert, C 4-snow/ice) C Iw1(Iv) : Iv=1 to Nintv; starting index for spectral summation C Iw2(Iv) : Iv=1 to Nintv; ending index for spectral summation C Jdprev : value of -Jday- in the previous entry C Misval : number representing missing data C Maxcld : maximum number of cloud layers C Maxout : maximum number of user-defined output levels C Nfpar : number of flux parameters C Nintv : no. of spectral intervals for output C Nlcld : number of cloud layers C Nout : number of user-defined output levels not including top C and surface C Nscene : no. of scene types C Nwl : no. of spectral intervals C Pole : true if polar region C Polnit : true if polar night, false otherwise C Pout(L) : L=1 to Nout; pressure (mb) of user-defined levels C not including top and surface C Prnt : if true, prints instantaneous values of fluxes C in file -Fnlst- C C Smumax : max value of solar zenith cosine to be considered C Smumin : min value of solar zenith cosine to be considered C Solcon : solar constant (W/m^2) C C Sundis : sun-earth distance factor C C Taucld(Is) : Is=1 to Nscene; optical depth C Vernum : version number of Modflx C Xdiff(it) : it=1, nintv; weight of diffuse surface albedo (0-1) C C +-------------------------------------------------------------------+ C Calling Tree (omitting calls to Errmsg): C C Modflx-+-Chkdim-+-Wrtdim C +-Chkinp-+-Wrtbad C +-Sundat C +-Getmod C +-Auxdat C +-Setcth C +-Solrad-+-Wrtdim C +-Prtins C +-------------------------------------------------------------------+ class(*) ,intent(in) :: aerObj !yma 160410 real :: Solcon C C .. Parameters .. INTEGER Maxcld, Maxlev, Maxout, Nwl, Maxint,Maxfpr, & Maxtyp, Maxwvl PARAMETER ( Nwl=7, Maxtyp=3, & Maxwvl=7,Maxfpr=5,Maxint=3, & Maxlev=50, Maxout=8, Maxcld=1 ) C .. C .. Scalar Arguments .. REAL Glat, Glon, Misval, Ozone, Pwater, & Solamu, Smumax, Smumin, Time, Tauaer, Cwcov, Cicov, & Tauwcd, Tauicd, Cldwer, Cldier, Pcldt, Psrf, Zsrf INTEGER Nwvl, Jday, Month LOGICAL Prnt CHARACTER(256) Header !CHARACTER (80) Header C .. Array Arguments .. REAL Flxall(Maxfpr,Maxint), Flxclr(Maxfpr,Maxint), & Albdir(Maxwvl), Albdif(Maxwvl) C .. C .. Local Scalars .. REAL :: Declin, Sundis, Pi INTEGER In, Ip, It, Isrtyp, Nout, Nzc, & Nfpar, Nintv, Nscene, Iaer, Iatm INTEGER Jdprev LOGICAL Badinp, Dimerr, Newday, Polnit,Climt LOGICAL First, delta CHARACTER (256) Fnls, Fnvt SAVE Jdprev, First C .. C .. Local Arrays .. REAL Fdntot( Maxlev ), Fdrtot( Maxlev ), & Fuptot( Maxlev ), Pc( Maxlev ), & Pout( Maxout ), Refdif( Maxwvl ), Refdir( Maxwvl ), & Zcldb( Maxcld ), Zcldt( Maxcld ) REAL Frac(Maxtyp),Clder(Maxtyp), & Tauclm(Maxtyp), Taucld(Maxtyp) REAL Flx(Maxfpr,Maxint,Maxtyp),xdiff(maxtyp) INTEGER Iw1(Maxint), Iw2(Maxint),Nlcld(Maxtyp) integer :: Jdayp real :: Vernum C .. Data statements .. DATA First / .TRUE. / , Vernum / 1.0 / , & Jdayp / -1 / SAVE pi, sundis, declin C .. Intrinsic Sunroutine .. INTRINSIC ASIN C .. External Subroutines .. c EXTERNAL Auxdat, Chkinp, Chkdim, Getmod, c & Prtins, Setcth, Solrad, Sundat Nintv=3 Nfpar=5 Nscene=3 Nout = 0 Pout = (/ 30., 180., 310., 440., 560., 680., 800., 1000. /) xdiff(1)=0.05 xdiff(2)=0.95 xdiff(3)=0.95 C .. IF ( First ) THEN Pi = 2. * Asin( 1.0 ) C ** Set up path and names of vegetation type and C ** cloud layer thickness model coefficients, and C ** output (listing) file C Fnvt = 'Auxdata/sfc/'//'vegtyp.dat' Fnls = 'Modflux.out' OPEN ( Unit=16, File=Fnls, Form='FORMATTED', Status='UNKNOWN' ) C ** Check dimensions and print un-changing inputs C CALL Chkdim( Smumin, Smumax, Solcon, Header, Vernum, Prnt, & Misval, Nwvl, Maxwvl,Iw1,Iw2 ) First = .FALSE. END IF !--------------------->>> userAer select type(aerObj) type is (real) Tauaer = aerObj type is (SRB_UserAer_Form1) Tauaer = aerObj%Tauaer type is (SRB_UserAer_Form2) Tauaer = aerObj%Tauaer end select !--------------------->>><< Return C IF ( Pwater .NE. Misval .AND. Ozone .NE. Misval ) RETURN C ** TROPICAL MODEL IF (Glat .GE. -25.0 .AND. Glat .LE. +25.0) THEN O3amt = 0.246 H2oamt = 4.12 GO TO 10 END IF C ** MIDLATITUDE MODEL IF (Abs(Glat) .GT.25.0 .AND. Abs(Glat) .LE.55.0) THEN IF (Glat .GT.0.0 .AND. Amjjas .OR. Glat. LT.0.0 & .AND. Jfmond) THEN C ** SUMMER O3amt = 0.318 H2oamt = 2.92 ELSE C ** WINTER O3amt = 0.396 H2oamt = 0.85 ENDIF GO TO 10 ENDIF C ** SUBARCTIC MODEL IF (Abs(Glat).GT.55.0) THEN IF (Glat.GT. 0.0 .AND. Amjjas .OR. Glat .LT. 0.0 & .AND. Jfmond) THEN C ** SUMMER O3amt = 0.344 H2oamt = 2.09 ELSE C ** WINTER O3amt = 0.478 H2oamt = 0.42 ENDIF ENDIF C C C ** REPLACE MISSING DATA FROM CLIMATOLOGY 10 IF (Pwater .LT. 0.0) Pwater = H2oamt IF (Ozone .LT. 0.0) Ozone = O3amt RETURN END C************************************************************************* SUBROUTINE Chkdim( Smumin, Smumax, Solcon, Header, Vernum, Prnt, & Misval, Nwlint, Maxwvl,Iw1,Iw2 ) C*********************************************************************** C C Checks dimensions and prints input variables that do not change C during a run C C INPUT VARIABLES: C C Header : 256-character long description C Maxwvl : max no. of spectral intervals in Subroutine Atsk13 C Misval : special value representing missing data C Nwlint : number of spectral intervals in Subroutine Atsk13 C Prnt : flag to turn on printing of input data and fluxes C Smumax : maximum solar zenith angle cosine to be used C Smumin : minimum solar zenith angle cosine to be used C Solcon : solar constant (W/m^2) for the total spectrum C Vernum : version number of Atsk13 C C OUTPUT VARIABLES: C Iw1(Iv) : Iv=1 to Nintv; starting index for spectral summation C Iw2(Iv) : Iv=1 to Nintv; ending index for spectral summation C C Called by- Atsk13 C Calls- Errmsg, Wrtdim C +-------------------------------------------------------------------+ C C C .. Scalar Arguments .. CHARACTER Header * ( * ) LOGICAL Prnt INTEGER Maxwvl, Nwlint REAL Misval, Smumax, Smumin, Solcon, Vernum INTEGER Iw1(*), Iw2(*) C .. C .. Local Scalars .. LOGICAL Dimerr C .. c LOGICAL Wrtdim C .. External Subroutines .. c EXTERNAL Errmsg, Wrtdim C .. C ** Set up starting and ending indexes for C ** output spectral intervals (1=sw, 2=visible, C ** 3=nir) Iw1(1) = 1 Iw2(1) = 7 Iw1(2) = 2 Iw2(2) = 4 Iw1(3) = 5 Iw2(3) = 7 C .. Dimerr = .FALSE. IF ( Maxwvl .LT. Nwlint ) Dimerr = Wrtdim( 'Maxwvl', Nwlint ) IF ( Dimerr ) CALL Errmsg( 'ATSK13--Dimension errors', .TRUE. ) WRITE ( 16, FMT=9000 ) Vernum WRITE ( 16, FMT='(/,2A)' ) ' (1) HEADER= ', Header WRITE ( 16, FMT=9010 ) Prnt WRITE ( 16, FMT=9020 ) Smumin, Smumax, Solcon, Misval RETURN 9000 FORMAT ( 1X, 125 ('*'), /, 34X, & 'MODIS SHORTWAVE RADIATION BUDGET ALGORITHM, VERSION ', & F5.2, /, 1X, 125 ('*') ) 9010 FORMAT ( /, ' (2) OPTIONS:', /, 5X, 'PRNT=', L2 ) 9020 FORMAT ( /, ' (3) PARAMETERS:', /, 5X, 'SMUMIN=', F5.2, $ ', SMUMAX=', F5.2, ', SOLCON=', F6.1, ', MISVAL=', & F8.1 ) END C************************************************************************** SUBROUTINE Chkinp( Smumin, Smumax, Solcon, Glat, Glon,Time, & Jday, Month,Cwcov,Cicov, Pwater, Ozone, & Solamu,Tauaer, Cldwer,Cldier,Pcldt, & Tauwcd,Tauicd, Albdir, Albdif, Misval, & Nwlint, Nwvl ) C C Checks input data C C INPUT VARIABLES: for a description of input variables see C Subroutine Atsk13 C C OUTPUT VARIABLES: none C C Called by- Atsk13 C Calls- Errmsg, Wrtbad C +-------------------------------------------------------------------+ C C .. Scalar Arguments .. INTEGER Jday, Nwlint, Nwvl,Month REAL Cwcov,Cicov, Glat, Glon, Misval, Ozone, Pwater, & Smumax, Smumin,Cldwer, Cldier, Pcldt, & Solamu, Solcon, Tauaer, Tauwcd, Tauicd, Time C .. C .. Array Arguments .. REAL Albdif( * ), Albdir( * ) C .. C .. Local Scalars .. LOGICAL First, Inperr INTEGER Iw C .. C c LOGICAL Wrtbad C .. External Subroutines .. c EXTERNAL Errmsg, Wrtbad C .. C .. Save statement .. SAVE First C .. C .. Data statements .. DATA First / .TRUE. / C .. Inperr = .FALSE. C ** In the first pass, check variables that do C ** not change during a run C IF ( First ) THEN First = .FALSE. IF ( Smumin .LT. 0.0 ) Inperr = Wrtbad( 'Smumin' ) IF ( Smumax .LT. Smumin .OR. Smumax .GT. & 1.0 ) Inperr = Wrtbad( 'Smumax' ) IF ( Solcon .LT. 0.0 ) Inperr = Wrtbad( 'Solcon' ) IF ( Nwvl .NE. Nwlint ) Inperr = Wrtbad( 'Nwvl' ) END IF C ** Check for valid range of coordinates, day and C ** time of day. C IF ( Glat .LT. -90. .OR. Glat .GT. 90. ) Inperr = Wrtbad( 'Glat' ) IF ( Glon .LT. -180. .OR. Glon .GT. 180.)Inperr = Wrtbad( 'Glon' ) IF ( (Time .LT. 0. .OR. Time .GT. 24.) .and. & Time /= 100.65 ) Inperr = Wrtbad( 'Time' )!130204,MaYingtao,time=9999 for daily 1deg processing. !IF ( Time .LT. 0. .OR. Time .GT. 24. ) Inperr = Wrtbad( 'Time' ) IF ( Jday .LT. 1 .OR. Jday .GT. 366 ) Inperr = Wrtbad( 'Jday' ) IF ( Month .LT. 1 .OR. Month .GT. 12 ) Inperr = Wrtbad( 'Month' ) C ** Check for valid values of satellite-observed C ** and derived parameters. A missing value of C ** cloud cover is assumed to indicate a missing C ** cell. All parameters must exist if the cloud C ** cover is not a missing value. Exceptions: C ** aerosol/cloud optical depth can be missing C ** when the scene is totally cloudy/clear. C IF ( Cwcov .NE. Misval .AND. (Cwcov .LT. 0.0 .OR. & Cwcov .GT. 1.001) ) Inperr = Wrtbad( 'Cwcov' ) IF ( Cicov .NE. Misval .AND. (Cicov .LT. 0.0 .OR. & Cicov .GT. 1.0) ) Inperr = Wrtbad( 'Cicov' ) IF (Cwcov .NE. Misval .AND. Cicov .NE. Misval .AND. & Cwcov+Cicov .GT. 1.001) Inperr = Wrtbad( 'Cover' ) DO 10 Iw = 1, Nwlint IF ( ( Cwcov .NE. Misval .OR. Cicov .NE. Misval) & .AND. (Albdir(Iw).EQ.Misval.OR. & Albdif(Iw).EQ.Misval.OR.Albdir(Iw).LT.0.0.OR. & Albdif(Iw).LT.0.0.OR.Albdir(Iw).GT.1.0.OR. & Albdif(Iw).GT.1.0) ) Inperr = Wrtbad( 'Albdir/Albdif' ) 10 CONTINUE IF ( (Cwcov .NE. Misval .OR. Cicov .NE. Misval) & .AND. (Solamu .EQ. Misval .OR. Solamu .LT. 0.0 .OR. & Solamu .GT. 1.0) ) Inperr = Wrtbad( 'Solamu' ) IF ( Cwcov .NE. Misval .AND. Cicov .NE. Misval .AND. & Cwcov+Cicov .NE. 1.0 .AND. & Tauaer .LT. 0.0 ) Inperr = Wrtbad( 'Tauaer' ) IF ( Cwcov .NE. Misval .AND. Cicov .EQ. Misval .AND. & Cwcov .NE. 1.0 .AND. Tauaer .LT. 0.0) & Inperr = Wrtbad( 'Tauaer' ) IF ( Cicov .NE. Misval .AND. Cwcov .EQ. Misval .AND. & Cicov .NE. 1.0 .AND. Tauaer .LT. 0.0) & Inperr = Wrtbad( 'Tauaer' ) IF ( Cwcov .NE. Misval .AND. Cwcov .NE. 0.0 .AND. & Tauwcd .LT. 0.0 ) Inperr = Wrtbad( 'Tauwcd' ) IF ( Cicov .NE. Misval .AND. Cicov .NE. 0.0 .AND. & Tauicd .LT. 0.0 ) Inperr = Wrtbad( 'Tauicd' ) IF ( Cwcov .NE. Misval .AND. Cwcov .NE. 0.0 .AND. & Cldwer .LT. 0.0 ) Inperr = Wrtbad( 'Cldwer' ) IF ( Cicov .NE. Misval .AND. Cicov .NE. 0.0 .AND. & Cldier .LT. 0.0 ) Inperr = Wrtbad( 'Cldier' ) IF ( ((Cwcov .NE. Misval .AND. Cwcov .NE. 0.0) .OR. & (Cicov .NE. Misval .AND. Cicov .NE. 0.0)) & .AND. Pcldt .LT. 0.0) Inperr = Wrtbad( 'Plcdt' ) c IF ( Pcldt .NE. Misval .AND. Pcldt .LT. c & 0.0 ) Inperr = Wrtbad( 'Pcldt' ) IF ( Pwater .NE. Misval .AND. Pwater .LT. & 0.0 ) Inperr = Wrtbad( 'Pwater' ) IF ( Ozone .NE. Misval .AND. Ozone .LT. & 0.0 ) Inperr = Wrtbad( 'Ozone' ) IF ( Inperr ) CALL Errmsg( 'ATSK13-Bad input data', .TRUE. ) RETURN END C********************************************************************** SUBROUTINE Setcth( Lat, Month, Isrftp,ctp, Pcldt, & Cldtau,Zsrf, Zcldt, Zcldb) C*********************************************************************** C Estimate cloud base height using a statistical model based on C 20 years climatology from rawinsonde observations(Wang et al., C 2000). Necessary inputs are cloud top height, month of year and C geolocations. C C INPUT: c Ctp: cloud type(1:water;2:ice) C Fnct: File containing cloud layer thickness model coefficient C Lat: Latitude C Month: Month of year (Northern Hemisphere) C Isrftp: surface type; 1=ocean; 2=land; 3=desert, 4=snow/ice C Srftp: Surface index; 1: Ocean, 2: Land C Pcldt: Cloud top height (Pressre Mpa) C Cldtau:Cloud optical thickness C Zsrf: Surface elevation C OUTPUT C Zcldt: Cloud top height( m ) C Zcldb: Cloud base height( m ) C C .. Implicit None Statement .. IMPLICIT NONE C .. C .. C .. Scalar Arguments .. LOGICAL First INTEGER Month, Srftp,Isrftp, ctp REAL Zcldt, Zcldb, Pcldt, Cldtau, Lat,Zsrf CHARACTER (256) Fnct C .. C .. Local Scalars .. INTEGER NZ,NPLBTB, Latidx, Il, Lindex,Ios REAL Ht, Cldthk, ztop, zdown C .. Local arrays REAL THKCOE(7,2,9,12),PLBTAB(12),ZLBTAB(12,72,12) C .. Save statement .. SAVE First, THKCOE, PLBTAB, ZLBTAB, NPLBTB C .. C .. Data statements .. DATA First / .TRUE. / C .. Fnct = 'Auxdata/'//'cld_thk.data' IF ( First ) THEN First = .FALSE. C ** Read coeffcients of cloud layer thickness model C OPEN(UNIT=8,FILE=Fnct, * STATUS='old',ACCESS='SEQUENTIAL',FORM='UNFORMATTED', * IOSTAT=IOS,convert='Big_Endian') Read(8,err=9)THKCOE,NPLBTB,PLBTAB,ZLBTAB CLOSE (8) GO TO 20 9 STOP 9 20 CONTINUE END IF C :: Obtain surface index Srftp=1 IF (isrftp .NE. 1) Srftp=2 C :: Obtain cloud top height by interpolating cloud top pressure into C cloud layer thickness model pressure levels. Latidx= INT((Lat+90.)/2.5)+1 IF (Latidx .GT. 72.0) Latidx=Latidx-1 C :: check cloud top pressure IF (ctp .EQ. 1)THEN IF (Pcldt .GE. 950. .OR. Pcldt .LT. 450.) THEN Pcldt=750. ENDIF ELSE IF (Pcldt .GE. 500.) THEN Pcldt=250. ENDIF ENDIF DO 10 Il=1,11 IF ((Pcldt .LT. PLBTAB(Il)) .AND. & (Pcldt .GE. PLBTAB(Il+1))) THEN Zdown=ZLBTAB(Il, Latidx, Month) Ztop=ZLBTAB(Il+1, Latidx, Month) Lindex=Il ENDIF 10 CONTINUE Zcldt=Zdown +(Ztop-Zdown) $ *(PLBTAB(Lindex)-Pcldt)/(PLBTAB(Lindex)-PLBTAB(Lindex+1)) C :: Obtain latutude zone index. If in south hemisphere, get month of year. Nz=INT(ABS(Lat)/10.0) IF (AMOD(ABS(Lat),10.0) .GT. 0.0) Nz=Nz+1 IF (Lat .LT. 0) Month=Month+6 IF (Month .GT. 12) Month=Month-12 C :: Determine if cloud belongs to deep convective clouds( cldtau >23 and C pcldt > 440 mpa). For deep convective cloud, put cloud base height C 800 m above surface IF ((Cldtau .GT. 23.0) .AND. (Pcldt .LT. 440.0))THEN Zcldb=Zsrf+800. ELSE cldthk = THKCOE(1,Srftp,Nz,Month) + $ Zcldt*THKCOE(2,Srftp,Nz,Month) $ + Zcldt**2*THKCOE(3,Srftp,Nz,Month) $ + Zcldt**3*THKCOE(4,Srftp,Nz,Month) Zcldb=Zcldt-cldthk C :: Check unreasonable estimate and replace it with mean value IF ((Zcldt .GT. 7000.) .AND. (cldthk .GT. Zcldt) $ .OR. (cldthk .LT. 0.0))THEN Zcldb=Zcldt-THKCOE(5,Srftp,Nz,Month) ENDIF IF (cldthk .LT. THKCOE(6,Srftp,Nz,Month))THEN Zcldb=Zcldt-THKCOE(7,Srftp,Nz,Month) ENDIF IF (Zcldb .LT. THKCOE(7,Srftp,Nz,Month))THEN Zcldb=THKCOE(7,Srftp,Nz,Month) ENDIF IF (Zcldb .LT. Zsrf) THEN Zcldb=Zsrf+400. ENDIF IF (Zcldb .GT. Zcldt) THEN Zcldt=Zcldb+1000. ENDIF ENDIF END C :: End of subroutine Setcth C**********************************************************************C SUBROUTINE Getmod( Glat, Glon, Flname, Iaer, Isrtyp) C**********************************************************************C C C Reads in vegetation type data of Matthews; associates surface type C with albedo models of Briegleb et al.; selects appropriate aerosol C model. C C References: 1) Briegleb, B. P., P. Minnis, V. Ramanathan and C E. Harrison, 1986: Comparison of regional clear sky C albedos inferred from satellite observations and C model calculations. J. Climate Appl. Meteor., C 25, 214-226. C C 2) Matthews, E., 1985: Atlas of archived vegetation, C land-use and seasonal albedo data sets. NASA C Technical Memorandum 86199, February 1985. C C C I N P U T: C C Flname : path and name of file containing surface type data C Glat : latitude of gridpoint (-90(S) to +90(N)) C Glon : longitude of gridpoint (-180(W) to +180(E)) C C O U T P U T : C C Iaer : aerosol model id (1=Cont; 2=Desert; 3=Mari; C 4=Arctic; 5=Antarctic, ) C Isrtyp : surface type; 1=ocean; 2=land; 3=desert, 4=snow/ice C Tauclm : climatological value of aerosol optical depth at C 0.55 microns C C I N T E R N A L V A R I A B L E S : C C First : true on first entry, false thereafter C Frslat : latitude of center of 1st cell of matthews' grid C Frslon : longitude of center of 1st cell of matthews' grid C Ialbmo(Iv) : Iv=1 to 32, array containing the albedo model id C of Briegleb et al. C Iaer : aerosol id associated with the surface type C (1=cont; 2=desert; 3=mar; 4=arctic;5=antarctic) C Iaerp : value of -Iaer- in previous pass C Ilat : latitude index of gridpoint in Matthews' maps C Ilon : longitude index of gridpoint in Matthews' maps C Isrmod : surface albedo model id: C -1 - snow/ice C 1 - mixed farming, tall grassland C 2 - tall/medium grassland, shrubland C 3 - short grassland, meadow and shrubland C 4 - evergreen forest C 5 - mixed deciduous forest C 6 - deciduous forest C 7 - tropical evergreen broadleaved forest C 8 - medium/tall grassland, woodlands C 9 - desert C 10 - tundra C 11 - water C Isrmp : value of -Isrmod- in previous pass C Newaer : true if current aerosol model differs from the C previous one, false otherwise C Newatm : true if current atmospheric model differs from the C previous one, false otherwise C Vegtyp(J,I) : J=1 to 180, I=1 to 360, vegetation type of an C 1 degree by 1 degree (lon/lat) area C C C .. Scalar Arguments .. CHARACTER Flname * ( * ) INTEGER Iaer, Ilat, Ilon, Isrmod, Isrtyp REAL Glat, Glon, Tauclm C .. C .. Local Scalars .. LOGICAL First, Newaer, Newmod, Newscn INTEGER I, Iaerp, Isrmp, Isrtp, Iv, J REAL Frslat, Frslon, Taup C .. C .. Local Arrays .. INTEGER Ialbmo( 0:32 ) INTEGER (KIND=2) Vegtyp( 360, 180 ) C .. C .. Intrinsic Functions .. INTRINSIC Nint C .. C .. Save statement .. SAVE First, Isrtp, Iaerp, Isrmp, Taup, Vegtyp C .. C .. Data statements .. DATA First / .TRUE. / DATA Frslat / -89.5 / , Frslon / -179.5 / DATA Ialbmo / 11, 7, 7, 7, 4, 4, 8, 4, 4, 6, 5, 6, 8, 8, 8, 8, 8, & 2, 2, 3, 10, 3, 10, 8, 2, 3, 2, 2, 3, 3, 9, -1, 1 / C .. C C IF ( First ) THEN First = .FALSE. C C ** Initialize certain variables C Iaerp = -9999 Isrmp = -9999 Isrtp = -9999 C ** Read surface-type data of Matthews C OPEN ( Unit=7, File=Flname, Form='FORMATTED', Status='OLD' ) DO 10 J = 1, 180 READ ( 7, FMT='(250I2,110I2)', END=20, & ERR=30 ) ( Vegtyp(I,J), I=1, 360 ) 10 CONTINUE CLOSE ( Unit=7 ) END IF IF ( Glon.GT.180. ) Glon = Glon - 360. C C ** Get lat/lon indices for location C Ilat = Nint( Glat-Frslat+1. ) IF ( Ilat.GT.180 ) Ilat = 180 Ilon = Nint( Glon-Frslon+1. ) IF ( Ilon.GT.360 ) Ilon = 360 C C ** Find surface type for location C Iv = Vegtyp( Ilon, Ilat ) C C ** Associate surface type with one of the C ** land albedo models of Briegleb et al. C Isrmod = Ialbmo( Iv ) Newmod = Isrmod .NE. Isrmp Isrmp = Isrmod IF ( Newmod ) THEN C ** New albedo model: Check if ocean, vegetative C ** land, desert or snow/ice scene C Isrtyp = 2 IF ( Isrmod.EQ.11 ) THEN Isrtyp = 1 ELSE IF ( Isrmod.EQ.9 ) THEN Isrtyp = 3 ELSE IF ( Isrmod.EQ.-1 ) THEN Isrtyp = 4 END IF Newscn = Isrtyp .NE. Isrtp Isrtp = Isrtyp IF ( Newscn ) THEN C C ** New scene type: Assign aerosol model and C ** climatological aerosol optical depth C ** pertinent to ocean, vegetative land, desert C ** and snow/ice scene IF ( Isrtyp.EQ.1 ) THEN Iaer = 3 Tauclm = 0.128 ELSE IF ( Isrtyp.EQ.2 ) THEN Iaer = 1 Tauclm = 0.23 ELSE IF ( Isrtyp.EQ.3 ) THEN Iaer = 2 Tauclm = 0.5 ELSE IF ( Isrtyp.EQ.4 ) THEN Iaer = 3 Tauclm = 0.05 END IF Newaer = Iaer .NE. Iaerp Taup = Tauclm Iaerp = Iaer ELSE IF ( .NOT.Newscn ) THEN C C ** Use previous aerosol model for this location C Iaer = Iaerp Tauclm = Taup END IF ELSE IF ( .NOT.Newmod ) THEN C C ** Use previous values for this location C Isrtyp = Isrtp Iaer = Iaerp Tauclm = Taup END IF RETURN 20 WRITE(16, FMT='(//,1X,A)')'GETMOD--PREMAT END OF VEG TYPE FILE' 30 WRITE(16, FMT='(//,1X,A)')'GETMOD--ERROR READING VEG TYPE DATA' STOP END C************************************************************************C SUBROUTINE Prtins( Jday, Glat, Glon, Time, Solamu, Ozone, & Pwater,Flx,Flxall, Frac, Misval,Tauaer, & Nfpar,Nintv,Nscene, Zsrf, Psrf,Pcldt, & Cldwer, Cldier, Tauwcd, Tauicd ) C C Prints instantaneous flux, optical depth . (Normally C called only for debugging.) C C INPUT : see Subroutine Modflx for a description C OUTPUT : none C C Called by- Modflx C +-------------------------------------------------------------------+ real :: Tauaer, Tauwcd INTEGER Maxfpr, Maxint,Maxtyp PARAMETER ( Maxint=3, Maxfpr=5,Maxtyp=3 ) C C .. Scalar Arguments .. INTEGER Jday REAL Glat, Glon, Tauwld, Tauicd, Cldwer, Cldier, & Misval, Ozone, Pwater, Zsrf, Psrf, & Solamu, Time ,Pcldt, Tauear C .. Array Arguments .. REAL Flxall(Maxfpr,*), Frac(*), Taucld(3), & Flx(Maxfpr,Maxint,Maxtyp) ,Clder(3) INTEGER Nfpar,Nintv,Nscene C INTEGER :: Iscene(*) C .. Local Scalars .. LOGICAL First, Newday INTEGER Jdayp REAL Ccover integer :: Ip,In,It, Jdprev C .. C .. Save statement .. SAVE First C .. C .. Data statements .. DATA First / .TRUE. / , Jdprev / -1 / C .. C ** PRINT INSTANTANEOUS FLUX AND OPTICAL DEPTH IF (First) THEN First = .FALSE. WRITE (16,Fmt=90000) WRITE (16,Fmt=90010) END IF Newday = Jday .NE. Jdprev IF ( Newday ) THEN Jdprev = Jday WRITE (16,Fmt=90020) Jday END IF WRITE (16,Fmt=90030) Taucld(1)=Tauaer Taucld(2)=Tauwcd Taucld(3)=Tauicd Clder(1)=0.0 Clder(2)=Cldwer Clder(3)=Cldier DO 20 It = 1, Nscene IF ( It .EQ. 1 ) THEN WRITE (16,Fmt=90040) Glat, Glon, Time, Solamu, Ozone, & Pwater, Zsrf, Psrf,Pcldt,Taucld(It), Clder(IT), & ((Flx(Ip,In,It),Ip=1,Nfpar),In=1,Nintv), It-1, & Frac(It) ELSE WRITE (16,Fmt=90060)Taucld(It), Clder(IT), & ((Flx(Ip,In,It),Ip=1,Nfpar),In=1,Nintv), It-1, & Frac(It) END IF 20 CONTINUE Ccover = Frac(2) + Frac(3) WRITE (16,Fmt=90080) ((Flxall(Ip,In),Ip=1,Nfpar),In=1,Nintv), & Ccover RETURN 90000 FORMAT (//'1',183('*')/55X,'FLUX ESTMATES'/1X, & 183('*')/) 90010 FORMAT (/T2,'PCT--CLOUD TOP PRESSURE(MB)',T46, & 'LAT---LATITUDE (DEGREES)',T99, & 'SRFDDN-DIFFUSE SURFACE DOWNWARD FLUX (W/SQ M)'/T2, & 'CCOVR-CLOUD COVER (%)',T46,'LON---LONGITUDE (DEGREES)',T99, & 'SRFUP-SURFACE UPWARD FLUX (W/SQ M)'/T2, & 'CSZA--COS( SOLAR ZENITH ANGLE )',T46, & 'CER--CLOUD DROPLET EFFECTIVE RADIUS(MICRON)',T99, & 'TAU---OPTICAL DEPTH'/T2,'FRAC--FRACTIONAL COVERAGE',T46, & 'O3----OZONE AMOUNT (ATM-CM)',T99, & 'TOADN-TOA DOWNWARD FLUX (W/SQ M)'/T2, & 'GMT---NOMINAL GMT (HOURS)',T46, & 'SCENE-0=CLR; 1=WATER CLD; 2=ICE CLD',T99, & 'TOAUP-TOA UPWARD FLUX (W/SQ M)'/T2, & 'H2O---WATER VAPOR AMOUNT (CM)',T46, & 'SRFDN-SURFACE DOWNWARD FLUX (W/SQ M)',T99, & 'ZSF-SURFACE ELEVATION (KM)',T2, & 'PSF-SURFACE PRESSURE (MB)' ) 90020 FORMAT (//1X,'|',78('-'),'>',3X,'DAY OF YEAR=',I4,3X,'<', & 78('-'),'|') 90030 FORMAT (/65X,2X,'<--------- Shortwave --------->',2X, & '<---------- Visible ---------->',2X, & '<------- Near Infrared -------->'/1X, & ' LAT LON UTC CSZA O3 H2O ZSF PSF', & ' PCT TAU CER',3(' TOADN TOAUP SRFDN SRFUP SRFDDN') & ,' SCENE FRAC') C90040 FORMAT (1X,F6.2,F8.2,F6.2,3F6.3,F7.2,15F7.1,I7,F6.3) 90040 FORMAT (1X,F6.2,F8.2,F6.2,F5.2,F5.2,F5.1,F5.1,2F6.0,2F6.1, & 15F7.0,I7,F6.3) C90060 FORMAT (39X,F7.2,15F7.1,I7,F6.3) 90060 FORMAT (53X,2F6.1,15F7.0,I7,F6.3) 90080 FORMAT (65X,15F7.1,' MIXED',F6.3) END C************************************************************************** C**********************************************************************C SUBROUTINE Solrad( Musun, Rsun, Solcon, Iatm, Iaer, & Iw1, Iw2, Xdif, Refdir, Refdif, Zsrf, & Psrf, Ozone, Pwater, Nlcld, Zcldt, & Zcldb, Taucld, Clder,Cldtp, aerObj, & Nout, Pout, Maxcld, Maxout, Fdrtot, & Fdntot,Fuptot, Flx, Nzc, Pc, Delta ) !yma 160410 C**********************************************************************C C ***With given atmosperic and surface parameters, using two-stream delta C ***Edington scheme to compute radiative flux at specified pressure level C C C C INPUT VARIABLES: C C Iaer : aerosol model id (1=Cont; 2=Desert; 3=Mari; C 4=Arctic; 5=Antarctic ) C Iatm : atmospheric model; one of six standard McClatchey C atmospheres: 1=Tropical, 2= Midlatitude summer, C 3=Midlatitude winter, 4=Subarctic summer, C 5=Subarctic winter, 6=US standard C Iw1(In) : In=1 to Nintv; starting index for spectral sum C Iw2(In) : In=1 to Nintv; ending index for spectral sum C Maxcld : maximum number of cloud layers C Maxout : maximum number of user-defined output levels C Musun : cosine of solar zenith angle C Nintv : no. of spectral intervals for output C Nlcld : number of cloud layers C Nout : number of user-defined output levels not including top C and surface C Ozone : total column amount of ozone (atm-cm) C Pout(L) : L=1 to Nout; pressure (mb) of user-defined levels not C including top and surface C Psrf : surface pressure (mb) C Pwater : user-defined total column amount of H2O (g/sq cm) C Refdif(I) : I=1 to 5; diffuse surface albedo in spectral interval -I- C Refdir(I) : I=1 to 5; direct surface albedo in spectral interval -I- C Rsun : correction factor for mean solar flux (square of ratio C of mean to actual earth-sun distance) C Solcon : solar constant (watts/m**2) C Taucld(L) : L=1 to Nlcld; optical depth of cloud layer -L- C Xdif : weight of diffuse surface albedo (0-1) C Zcldb(L) : L=1 to Nlcld; cloud-base height for cloud layer -L- (km) C Zcldt(L) : L=1 to Nlcld; cloud-top height for cloud layer -L- (hpa) C Zsrf : surface height (km) C C OUTPUT VARIABLES: C Flx C C INTERNAL VARIABLES: C C Cutaer : wavelength (microns) above which aerosol scattering is C to be neglected. ( default: 4.0 ) C Cutcld : wavelength (microns) above which cloud scattering is C to be neglected. ( default: 4.0 ) C Cutray : wavelength (microns) above which Rayleigh scattering is C to be neglected. ( default: 4.0 ) C Flux(Ip,In) :Ip=1,5;In=1,Nwl: spetral fluxes at surface and at top of atmosphere C Lambda : mean wavelength (microns) for current spectral C interval ( average of *Wl(Iwl)* and *Wl(Iwl+1)* ) C Noray : true, neglect Rayleigh scattering C Laer(Lc) : Lc=1 to Nzc, model-aerosol layer to be used in C computational layer Lc C Lcld(Lc) : Lc=1 to Nzc, model-cloud layer to be used in C computational layer Lc C Maxint : max no. of output spectral intervals C C Maxlev : max no. of vertical levels C Maxwvl : max no. of spectral intervals C Maxnum : max no. of terms in correlated-k distribution of C absorption C Nzc : number of vertical levels C Nwl : number of spectral intervals C Wavlen(I) : I=1 to Nwl+1; wavelength (micron); -Wavlen(I)- is the C lower boundary and -Wavlen(I+1)- is the upper boundary C of a spectral interval C +---------------------------------------------------+ C Calling tree C C Solrad-+-Setatm C +-Setaer C +-Setcld C +-Effamt C +-Scat C +-Absmol C +-Solar C +-Main2 C +-Abso3 C +-Prnflx C +----------------------------------------------------+ class(*) ,intent(in) :: aerObj !yma 160410 integer :: Maxcld, Maxout C .. Parameters .. INTEGER Maxwvl, Maxnum, Maxpar, Maxint,Maxlev PARAMETER ( Maxwvl=7, Maxnum=10, Maxint=3, Maxpar=5, & Maxlev=50 ) C .. C .. Scalar Arguments .. INTEGER Iaer, Ialb, Iatm, Isnowc, Isrf, Maxcld & Maxout, Nlcld, Nout, Nzc REAL Musun, Ozone, Psrf, Rsun, Solcon, Tauaer, Xalb, Xdif, & Zsrf, Pwater C .. C .. Array Arguments .. INTEGER Cldtp(Maxcld), Iw1(Maxint), Iw2(Maxint) REAL Fdntot( Maxlev ), Fdrtot( Maxlev ), & Fuptot( Maxlev ), & Pout( Maxout ), Refdif( Maxwvl ), Refdir( Maxwvl ), & Taucld( Maxcld ), Clder(Maxcld), & Pcldt(Maxcld),Zcldb( Maxcld ), Zcldt( Maxcld ), & Flx(Maxpar, Maxint),Pc( Maxlev ) C .. C .. Local Scalars .. LOGICAL Noaer, Nocld, Nomabs, Noray, Prnspr, Prnt INTEGER Iwl, Lc, Nab, Nkn, Nuvis, Nwl, Ip, Il, Nintv REAL Alb, Albdif, Albdir, Cutaer, Cutcld, Cutray, Flxsun, Lambda, & Solflx, Xh2o C .. C .. Local Arrays .. INTEGER Laer( Maxlev-1 ), Lcld( Maxlev-1 ), Lprn( Maxlev ) REAL Co3( Maxlev ), Dtauex( Maxlev-1 ), Dtausc( Maxlev-1 ), & Dteaer( Maxlev-1, Maxwvl ), Dtecld( Maxlev-1, Maxwvl ), & Fld( Maxlev ), Fldn( Maxlev ), Flup( Maxlev ), & Gf( Maxlev-1 ), Gfaer( Maxlev-1, Maxwvl ), & Gfcld( Maxlev-1, Maxwvl ), H2oden( Maxlev ), & H2oeff( Maxlev-1 ), Kn( Maxnum ), & Omaer( Maxlev-1, Maxwvl ), & Omcld( Maxlev-1, Maxwvl ), Pkn( Maxnum ), Qv( Maxlev ), & T( Maxlev ), Wavlen( Maxwvl+1 ), & Xl( Maxlev ), Xlstar( Maxlev ), Zc( Maxlev ), & Flux(Maxwvl,Maxpar) LOGICAL delta C .. C .. External Subroutines .. c EXTERNAL Absmol, Abso3, Effamt, Main2, Prnflx, Scat, c & Setaer, Setatm, Setcld, Solar,Sums C .. C .. Data statements .. DATA Nwl / 7 / , Nuvis / 4 /, Nintv / 3 / DATA Wavlen / 0.2, 0.4, 0.5, 0.6, 0.7, 1.19, 2.38, 4.0 / DATA Cutray / 4.0 / , Cutaer / 4.0 / , Cutcld / 4.0 / DATA Prnspr / .False. / , Prnt / .FALSE. / !140102 DATA Nomabs / .FALSE. / C .. C ** Set up atmospheric vertical profile C CALL Setatm( Prnt, Iatm, Iaer, Maxlev, Ozone, Pwater, & Zsrf/1000., Psrf, Nout, Pout, Nlcld, & Zcldt/1000., Zcldb/1000., Nzc, Zc, & Pc, T, Qv, H2oden, Co3, Xh2o, Laer, Lcld, Lprn, & aerObj ) ! 160410 C ** Compute aerosol extinction optical depth, single C ** scattering albedo, and asymmetry factor C CALL Setaer( Nzc, Zc, Iaer, aerObj, Laer, Nwl, Maxlev, Maxwvl, & Dteaer, Omaer, Gfaer ) ! 160410 C ** Compute cloud extinction optical depth, single C ** scattering albedo, and asymmetry factor C CALL Setcld( Nzc, Zc, Nlcld, Cldtp,Taucld,clder, Zcldt, & Zcldb, Lcld, Nwl, & Wavlen, Nuvis, Musun, Maxlev, Maxwvl, Dtecld, Omcld, & Gfcld ) C ** Compute effective absorber amounts of water vapor C ** and ozone C CALL Effamt( Nzc, Zc, Pc, T, H2oden, Xh2o, Co3, Musun, Maxlev, & H2oeff, Xl, Xlstar ) !---------------->>> userAer select type(aerObj) type is (real) Tauaer = aerObj type is (SRB_UserAer_Form1) Tauaer = aerObj%Tauaer type is (SRB_UserAer_Form2) Tauaer = aerObj%Tauaer end select !---------------->>><<< userAer Fdrtot = 0.0 Fdntot = 0.0 Fuptot = 0.0 DO 20 Iwl = 1, Nwl Lambda = 0.5 * ( Wavlen(Iwl)+Wavlen(Iwl+1) ) Noray = Lambda .GT. Cutray Noaer = Lambda .GT. Cutaer .OR. Tauaer .EQ. 0.0 Nocld = Lambda .GT. Cutcld .OR. Nlcld .EQ. 0 C ** For combined Rayleigh and particulate scattering, C ** calculate extinction and scattering optical depths, C ** and asymmetry factors CALL Scat( Iwl, Nzc, Zc, Pc, Qv, Noray, Noaer, Nocld, Maxlev, & Dteaer(1,Iwl), Omaer(1,Iwl), Gfaer(1,Iwl), & Dtecld(1,Iwl), Omcld(1,Iwl), Gfcld(1,Iwl), Dtauex, & Dtausc, Gf ) C C ** Obtain pseudo-molecular absorption coefficients C Nab = 0 IF ( .NOT.Nomabs ) CALL Absmol( Iwl, Maxnum, Nab, Nkn, Pkn, & Kn ) C ** User supplied albedo C Albdir = Refdir( Iwl ) Albdif = Refdif( Iwl ) C ** Weight direct and diffuse components of surface C ** albedo Alb = ( 1.-Xdif ) * Albdir + Xdif * Albdif C ** Calculate TOA solar irradiance C CALL Solar( Iwl, Musun, Rsun, Solcon, Solflx, Flxsun ) c C ** Solve for fluxes in this spectral interval as a C ** weighted sum of monochromatic multiple scattering C ** problems c CALL Main2( Nab, Nkn, Pkn, Kn, H2oeff, Nzc, Dtauex, Dtausc, Gf, & Maxlev, Musun, Alb, Solflx, Flxsun, Fld, Fldn, & Flup ,Delta) IF (.NOT. Delta)GOTO 200 C ** Account for ozone absorption C CALL Abso3( Lambda, Musun, Xl, Xlstar, Nzc, Maxlev, Rsun, & Solflx, Solcon, Fld, Fldn, Flup ) C ** Print fluxes for current spectral interval C IF ( Prnspr ) CALL Prnflx( Fld, Fldn, Flup, Musun, Nzc, Pc, & Maxlev, Wavlen(Iwl), Wavlen(Iwl+1) ) C ** Increment spectral sums of fluxes C DO 10 Lc = 1, Nzc Fdrtot( Lc ) = Fdrtot( Lc ) + Fld( Lc ) Fdntot( Lc ) = Fdntot( Lc ) + Fldn( Lc ) Fuptot( Lc ) = Fuptot( Lc ) + Flup( Lc ) 10 CONTINUE C ** Assign spectral fluxes at surface and at top of atmosphere Flux(Iwl,1)=Fldn(1) Flux(Iwl,2)=Flup(1) Flux(Iwl,3)=Fldn(Nzc) Flux(Iwl,4)=Flup(Nzc) Flux(Iwl,5)=Fldn(Nzc)-Fld(Nzc) c 20 CONTINUE C ** compute fluxes at output spectral intervals DO 40 Ip=1, Maxpar DO 30 Il=1, Nintv Flx(Ip,Il) = Sums( Flux(1,ip), Iw1(Il), Iw2(Il), Nwl) 30 CONTINUE 40 CONTINUE C ** Print total fluxes C IF ( Prnt ) CALL Prnflx( Fdrtot, Fdntot, Fuptot, Musun, Nzc, Pc, & Maxlev, Wavlen(1), Wavlen(Nwl+1) ) 200 RETURN END C*********************************************************************** SUBROUTINE Absmol( Iwl, Maxnum, Nab, Nkn, Pkn, Kn ) C*********************************************************************** C C For current spectral interval, finds pseudo-molecular C absorption coefficients C Modified Lacis and Hansen H2O parameterization C k-distribution parameterization of H2O absorption, Pkn C and Kn( cm**2/g). (Ramaswamy, and S. M. Freidenreich, 1992: C A study of broadband parameterizations of the solar radiative C interactions with water vapor and water drops, J. Geophys. C Res., D11, 11,487-11,512) C .. Implicit None Statement .. IMPLICIT NONE C .. C .. Parameters .. INTEGER Nkntab, Nwltab PARAMETER ( Nkntab=10, Nwltab=7 ) C .. C .. Scalar Arguments .. INTEGER Iwl, Maxnum, Nab, Nkn C .. C .. Array Arguments .. REAL Kn( Maxnum ), Pkn( Maxnum ) C .. C .. Local Scalars .. LOGICAL Dimerr INTEGER K,L C .. C .. Local Arrays .. INTEGER Numabs( Nwltab ) REAL Kntab( Nkntab ,2), Pkntab( Nkntab ,4) C .. C .. External Functions .. c LOGICAL Wrtdim c EXTERNAL Wrtdim C .. C .. External Subroutines .. c EXTERNAL Errmsg C .. C .. Intrinsic Functions .. INTRINSIC Sum C .. C .. Data statements .. DATA (Pkntab(L,1), L=1,10) / 7.332E-1, 2.199E-1, 2.4611E-2, & 1.3891E-2, 6.90802E-3, 7.96458E-4, 2.08745E-4, & 1.75948E-4, 1.57633E-4, 8.54838E-5 / DATA (Pkntab(L,2), L=1,10) / 6.02392E-1, 1.78305E-1, 6.5137E-2, & 7.5077E-2, 4.37527E-2, 1.81407E-2, 7.68065E-3, & 5.0843E-3, 3.14907E-3, 1.28161E-3 / DATA (Pkntab(L,3), L=1,10) / 4.1872E-1, 1.18546E-1, 4.80756E-2, & 1.03762E-1,6.76036E-2, 8.32642E-2, 1.21417E-1, & 1.60241E-2, 1.70456E-2, 5.54177E-3 / DATA (Pkntab(L,4), L=1,10) / 1.00184E-1, 1.58381E-1, 1.306E-1, & 1.49868E-1,1.20244E-1, 6.57255E-2, 7.33715E-2, & 6.92753E-2, 1.13355E-1, 1.89953E-2 / DATA (Kntab(L,1), L=1,10) / .0000, .001, .0133, .0422, 0.1334, & 0.4217, 1.334, 5.623, 31.62, 177.8 / DATA (Kntab(L,2), L=1,10) / .001, .0133, .0422, 0.1334, 0.4217, & 1.334, 5.623, 31.62, 177.8, 1000. / DATA Numabs / 0, 0, 0, 1, 1, 1, 1/ C .. Dimerr = .FALSE. ! 160424 IF ( Nkntab.GT.Maxnum ) Dimerr = Wrtdim( 'Maxnum', Nkntab ) IF ( Dimerr ) CALL Errmsg( 'Absmol--Dimension error', .TRUE. ) Nab = Numabs( Iwl ) IF ( Nab.EQ.0 ) THEN Nkn = 0 Pkn = 0.0 Kn = 0.0 RETURN ELSE Nkn = Nkntab DO 10 K = 1, Nkntab Pkn( K ) = Pkntab( K, Iwl-3 ) IF ( Iwl.EQ.4 ) THEN Kn( K ) = Kntab( K, 1 ) ELSE Kn( K ) = Kntab( K, 2 ) ENDIF 10 CONTINUE END IF END C*********************************************************************** SUBROUTINE Abso3( Lambda, Musun, Xl, Xlstar, Nz, Maxlev, Rsun, & Solflx, Solcon, Fld, Fldn, Flup ) C*********************************************************************** C C Calculates ozone transmission in current spectral interval, C and attenuates flux for ozone absorption C C INPUT: C C Fld(L) : L=1 to Nz; downward direct flux at level -L- C Fldn(L) : L=1 to Nz; downward flux (diffuse+direct) at level -L- C Flup(L) : L=1 to Nz; same as *Fldn*, only upward C Lambda : wavelength (microns) C Musun : cosine of solar zenith angle C Nz : number of levels C Maxlev : maximum number of levels in calling program C Rsun : correction factor for mean solar flux (square of ratio C of mean to actual earth-sun distance) C Solcon : solar constant (watts/m**2) C Solflx : solar flux normal to beam (watts/m**2) C Xl(L) : L=1 to Nz; effective absorber path for downward C radiation (atm-cm) C Xlstar(L) : L=1 to Nz; effective absorber path for upward radiation C (atm-cm) C C OUTPUT: C C Fld(L) : L=1 to Nz; downward direct flux at level -L- C Fldn(L) : L=1 to Nz; downward flux (diffuse+direct) at level -L- C Flup(L) : L=1 to Nz; same as *Fldn*, only upward C LOCAL VARIABLES: C C Trdn(L) : L=1 to Nz; transmissivity between TOA and level -L- for C downward radiation C Trup(L) : L=1 to Nz; transmissivity between TOA and level -L- for C upward radiation C .. Scalar Arguments .. INTEGER Maxlev, Nz REAL Lambda, Musun, Rsun, Solcon, Solflx C .. C .. Array Arguments .. REAL Fld( Maxlev ), Fldn( Maxlev ), Flup( Maxlev ), Xl( Maxlev ), & Xlstar( Maxlev ) C .. C .. Local Scalars .. INTEGER L REAL Absdn, Absup, Frac C .. C .. Local Arrays .. REAL Trdn( Maxlev ), Trup( Maxlev ) C .. C .. External Functions .. c REAL Absuvi, Absvis c EXTERNAL Absuvi, Absvis C .. IF ( Lambda.LT.0.4 ) THEN DO 10 L = 1, Nz Absdn = Rsun * Solcon * Absuvi( Xl(L) ) Absup = Rsun * Solcon * Absuvi( Xlstar(L) ) Trdn( L ) = 1.0 - Absdn / Solflx Trup( L ) = 1.0 - Absup / Solflx 10 CONTINUE ELSE IF ( Lambda.GE.0.5 .AND. Lambda.LT.0.6 ) THEN Frac = 7.5E-3 * Musun + 0.6208 DO 20 L = 1, Nz Absdn = Frac * Rsun * Solcon * Absvis( Xl(L) ) Absup = Frac * Rsun * Solcon * Absvis( Xlstar(L) ) Trdn( L ) = 1.0 - Absdn / Solflx Trup( L ) = 1.0 - Absup / Solflx 20 CONTINUE ELSE IF ( Lambda.GE.0.6 .AND. Lambda.LT.0.7 ) THEN Frac = 7.5E-3 * Musun + 0.6208 DO 30 L = 1, Nz Absdn = ( 1.-Frac ) * Rsun * Solcon * Absvis( Xl(L) ) Absup = ( 1.-Frac ) * Rsun * Solcon * Absvis( Xlstar(L) ) Trdn( L ) = 1.0 - Absdn / Solflx Trup( L ) = 1.0 - Absup / Solflx 30 CONTINUE ELSE Trdn = 1.0 Trup = 1.0 END IF DO 40 L = 1, Nz Fld( L ) = Trdn( L ) * Fld( L ) Fldn( L ) = Trdn( L ) * Fldn( L ) Flup( L ) = Trup( L ) * Flup( L ) 40 CONTINUE END c*********************************************************************** REAL FUNCTION Absuvi( O3path ) c*********************************************************************** C C Computes the absorption of ozone in the ultraviolet (Hartley- C Huggins) band C C .. Scalar Arguments .. REAL O3path C .. Absuvi = ( 1.0820 *O3path/ (1.+138.6 *O3path)**0.805+ & 0.0658 *O3path/ (1.+ (103.6 *O3path)**3.0) ) END c*********************************************************************** REAL FUNCTION Absvis( O3path ) c*********************************************************************** C C Computes the absorption of ozone in the visible (Chappuis) band C C .. Scalar Arguments .. REAL O3path C .. Absvis = 0.02118 * O3path / ( 1.+.042 *O3path+ & .000323 *O3path**2. ) END c*********************************************************************** SUBROUTINE Delted( Solflx, Dttot, Omega, Gf, Maxlev, Alb, Rmu, & Nz, Dirflx, Difflx, Upwflx, Delta ) c*********************************************************************** C C Calculates up- and down-fluxes of radiation in a vertically C inhomogeneous atmosphere using the delta-Eddington C approximation ( this is a variant of the code given in C NCAR Technical Note TN-121+STR : 'The delta-Eddington C approximation for a vertically inhomogeneous atmosphere', C by W.J. Wiscombe (July, 1977) ) C C ** Note ** Requires the subroutine *Leqt1b* and *Leqt2b* from C the I.M.S.L. library C INPUT: C C Alb : surface albedo C Dttot(L) : L=1 to Nz; continuum extinction optical depth for layer L C Gf(L) : L=1 to Nz; asymmetry factor for layer L C Maxlev : max no. of levels in calling program C Nz : number of levels (level 1 = top, level nz =surface) C Omega(L) : L=1 to Nz; single scattering albedo for layer L C Solflx : incident beam flux at level 1 ( normal to beam ) C Rmu : cosine of solar zenith angle C C C OUTPUT: C C Dirflx(L) : L=1 to Nz; direct-beam flux at level -L- C ( includes some near-forward-scattered radiation C along with the actual direct beam ) C Difflx(L) : L=1 to Nz; diffuse down-flux at level -L- C ( this is less than the actual diffuse down-flux by C the same amount that the direct flux -Dirflx- is C augmented. ) C Upflx(L) : L=1 to Nz; diffuse up-flux at level -l- C C C LOCAL VARIABLES: C C Dir(L,Np) : L=1 to Nz, Np=1 to Nsol; direct-beam flux at level -L- C for Np-th sun angle C ( note -- because of the way the forward scattering C peak is truncated, *Dir* includes some C near-forward-scattered radiation along with C the actual direct beam ) C Flxd(L,Np): L=1 to Nz, Np=1 to Nsol; diffuse down-flux at level -L- C for Np-th sun angle C ( note -- this will be less than the actual diffuse C down-flux by the same amount that the C direct flux -Dir- is augmented. ) C Flxu(L,Np): L=1 to Nz, Np=1 to Nsol; diffuse up-flux at level -l- C for np-th sun angle C Prec : number to subtract from single-scattering albedo when it C is unity to avoid doing special case branching in the C method (ten times the machine precision) integer :: Maxlev integer, PARAMETER :: Nlde=50, Nlde2=2*Nlde-2, Maxsun=1 C .. Scalar Arguments .. INTEGER Nz REAL Alb, Rmu, Solflx C .. C .. Array Arguments .. REAL Difflx( Maxlev ), Dirflx( Maxlev ), Upwflx( Maxlev ) DOUBLE PRECISION Dttot( Maxlev-1 ), Gf( Maxlev-1 ), & Omega( Maxlev-1 ) C .. C .. Local Scalars .. LOGICAL Dimerr, delta INTEGER I, Ib, Ier, Ijob, Ip1, J, L, Last, Lastm2, Np, & Nsol, Nzm1 REAL A11, A12, A21, A22, C1, Cutpt, Det, Prec, Rmu0, & Solexp, T1, T2, T3, T4, T5, Xx DOUBLE PRECISION Ff C .. C .. Local Arrays .. REAL Dir( Nlde, Maxsun ),Flxd( Nlde, Maxsun ), Flxsun( Maxsun ), & Flxu( Nlde, Maxsun ), Musun( Maxsun ) DOUBLE PRECISION Alph( Nlde ), Beta( Nlde ), & Ex( Nlde ), & Exsun( Nlde ), & Lm( Nlde ), Pp( Nlde ), & Ss( Nlde2, 5 ), Tau( Nlde ), Tx( Nlde ), & Ty( Nlde ), Tz( Nlde ), Work( Nlde2, 10 ), & X( Nlde2 ) C .. C .. External Subroutines .. c EXTERNAL Leqt2b C .. C .. Intrinsic Functions .. INTRINSIC Abs, Epsilon, Exp, Sqrt C .. C .. Data statements .. DATA Cutpt / 10. / C .. C Prec = 10. * Epsilon( 1.0 ) delta=.true. Nsol = 1 Dimerr = .FALSE. IF ( Nsol.GT.Maxsun ) Dimerr = Wrtdim( 'Maxsun', Nsol ) IF ( Dimerr ) CALL Errmsg( 'Delted--Dimension error', .TRUE. ) Musun( 1 ) = Rmu C1 = 2. / 3. Nzm1 = Nz - 1 DO 10 I = 1, Nzm1 Lm( I ) = 0. Pp( I ) = 0. Ex( I ) = 0. Tx( I ) = 0. Ty( I ) = 0. Tz( I ) = 0. Tau( I ) = 0. Exsun( I ) = 0. Alph( I ) = 0. Beta( I ) = 0. 10 CONTINUE Tau( Nz ) = 0. C C ** Scale optical depth, sing-scat albedo, and C ** asymmetry factor and calculate various fcns of C ** these variables C Tau( 1 ) = 0. DO 20 I = 1, Nzm1 Ff = Gf( I )**2 Gf( I ) = Gf( I ) / ( 1.+Gf(I) ) Dttot( I ) = ( 1.-Omega(I) *Ff ) * Dttot( I ) Omega( I ) = ( 1.-Ff ) * Omega( I ) / ( 1.-Ff*Omega(I) ) IF ( Omega(I).GE.1.0 ) Omega( I ) = 1.0 - Prec T1 = 1. - Omega( I ) * Gf( I ) Lm( I ) = Sqrt( 3. * (1.-Omega(I)) *T1 ) Pp( I ) = C1 * Lm( I ) / T1 C ** Treat highly-absorbing layers, which would C ** otherwise cause ill-conditioning in the C ** penta-diagonal matrix "Ss", by assuming C ** the bottom boundary of the calculation is at C ** the top of the first such layer encountered IF ( Lm(I) *Dttot(I).GE.Cutpt ) Dttot( I ) = Cutpt / Lm( I ) Ex( I ) = Exp( Lm(I) *Dttot(I) ) T1 = Gf( I ) * ( 1.-Omega(I) ) Tx( I ) = 0.75 * Solflx * Omega( I ) * ( 1.+T1 ) Ty( I ) = 0.5 * Solflx * Omega( I ) Tz( I ) = 3. * T1 Tau( I+1 ) = Tau( I ) + Dttot( I ) 20 CONTINUE Ib = Nz C ** Set incident flux at top of atmosphere C DO 40 Np = 1, Nsol Flxsun( Np ) = Musun( Np ) * Solflx 40 CONTINUE IF ( Ib.LT.2 ) THEN DO 60 Np = 1, Nsol Dir( 1, Np ) = Flxsun( Np ) Flxd( 1, Np ) = 0. Flxu( 1, Np ) = 0. DO 50 I = 2, Nz Dir( I, Np ) = 0. Flxd( I, Np ) = 0. Flxu( I, Np ) = 0. 50 CONTINUE 60 CONTINUE ELSE IF ( Ib.EQ.2 ) THEN C C ** Do one-layer case separately, since multi-layer C ** formalism breaks down for one layer (system is C ** no longer penta-diagonal) C IF ( Abs(1.-Omega(1)).GT.Prec ) THEN C ** Absorbing case A11 = 1. - Pp( 1 ) A12 = 1. + Pp( 1 ) A21 = ( A12-Alb*A11 ) * Ex( 1 ) A22 = ( A11-Alb*A12 ) / Ex( 1 ) Det = A12 * A21 - A11 * A22 DO 70 Np = 1, Nsol Rmu0 = 1. / Musun( Np ) T1 = Rmu0**2 - Lm( 1 )**2 Alph( 1 ) = Tx( 1 ) / T1 Beta( 1 ) = Ty( 1 ) * ( Musun(Np) *Tz(1)+Rmu0 ) / T1 Solexp = Exp( -Rmu0*Tau(2) ) X( 1 ) = ( A12* (Alph(1)-Beta(1)+Alb* (Flxsun(Np)- & Alph(1)-Beta(1))) *Solexp- & A22* (Alph(1)+Beta(1)) ) / Det X( 2 ) = ( Alph(1)+Beta(1)-A11*X(1) ) / A12 Dir( 1, Np ) = Flxsun( Np ) Dir( 2, Np ) = Flxsun( Np ) * Solexp Flxd( 1, Np ) = A11 * X( 1 ) + A12 * X( 2 ) - Alph( 1 ) - & Beta( 1 ) Flxd( 2, Np ) = A11 * X( 1 ) * Ex( 1 ) + & A12 * X( 2 ) / Ex( 1 ) - & ( Alph(1)+Beta(1) ) * Solexp Flxu( 1, Np ) = A12 * X( 1 ) + A11 * X( 2 ) - Alph( 1 ) + & Beta( 1 ) Flxu( 2, Np ) = A12 * X( 1 ) * Ex( 1 ) + & A11 * X( 2 ) / Ex( 1 ) - & ( Alph(1)-Beta(1) ) * Solexp 70 CONTINUE ELSE C ** No absorption case (must be treated separately C ** to avoid overflows and also because fcnal forms C ** are quite different) T1 = ( 1.-Gf(1) ) * Tau( 2 ) T2 = 4. / 3. + ( 1.-Alb ) * T1 DO 80 Np = 1, Nsol T4 = 1. + 1.5 * Musun( Np ) T5 = 1. - 1.5 * Musun( Np ) Solexp = Exp( -Tau(2)/Musun(Np) ) X( 1 ) = ( 1.-Alb ) * ( T4+T5*Solexp ) / T2 X( 2 ) = T4 - C1 * X( 1 ) Dir( 1, Np ) = Flxsun( Np ) Dir( 2, Np ) = Flxsun( Np ) * Solexp T3 = 0.5 * Flxsun( Np ) Flxd( 1, Np ) = ( C1*X(1)+X(2)-T4 ) * T3 Flxd( 2, Np ) = ( (C1-T1) *X(1)+X(2)-T4*Solexp ) * T3 Flxu( 1, Np ) = ( -C1*X(1)+X(2)+T5 ) * T3 Flxu( 2, Np ) = ( - (C1+T1) *X(1)+X(2)+T5*Solexp ) * T3 80 CONTINUE END IF ELSE IF ( Ib.GT.2 ) THEN C C ** Multi-layer delta-eddington C ** The top and bottom boundary conditions plus C ** the flux continuity conditions at each C ** interior level form a penta-diagonal system C ** of 2*Nz-2 equations for the unknown constants C ** (2 for each layer). The columns of the "Ss" C ** array contain the diagonals of the coefficient C ** matrix, the lowermost diagonal in column 1, C ** etc.(This is the so-called band storage mode C ** required by IMSL routine Leqt2b). C Last = 2 * Ib - 2 Ss( 1, 1 ) = 0. Ss( 1, 2 ) = 0. Ss( 1, 5 ) = 0. Ss( Last, 1 ) = 0. Ss( Last, 4 ) = 0. Ss( Last, 5 ) = 0. Ss( 1, 3 ) = ( 1.-Pp(1) ) / Ex( 1 ) Ss( 1, 4 ) = ( 1.+Pp(1) ) * Ex( 1 ) Lastm2 = Last - 2 DO 90 J = 2, Lastm2, 2 I = J / 2 Ip1 = I + 1 Ss( J, 1 ) = 0. Ss( J, 2 ) = 1.0 Ss( J, 3 ) = 1.0 Ss( J, 4 ) = -1.0 / Ex( Ip1 ) Ss( J, 5 ) = -Ex( Ip1 ) Ss( J+1, 1 ) = -Pp( I ) Ss( J+1, 2 ) = Pp( I ) Ss( J+1, 3 ) = Pp( Ip1 ) / Ex( Ip1 ) Ss( J+1, 4 ) = -Pp( Ip1 ) * Ex( Ip1 ) Ss( J+1, 5 ) = 0. 90 CONTINUE Ss( Last, 2 ) = 1. + Pp( Ip1 ) - Alb * ( 1.-Pp(Ip1) ) Ss( Last, 3 ) = 1. - Pp( Ip1 ) - Alb * ( 1.+Pp(Ip1) ) C C ** Calculate the L-U decomposition of penta- C ** diagonal matrix "Ss" and save in "Work" C ** array C Ijob = 1 CALL Leqt2b( Ss, Last, 2, 2, Nlde2, X, 1, Nlde2, Ijob, Work, & Nlde2, Work(1,8), Ier ) IF (Ier .NE. 0) GO TO 150 C C ** For each sun angle, calculate the r.h.s.of C ** the banded system, solve, and use the C ** solution to construct the fluxes at each C ** level C Ijob = 2 DO 120 Np = 1, Nsol Rmu0 = 1. / Musun( Np ) T1 = Rmu0**2 - Lm( 1 )**2 Alph( 1 ) = Tx( 1 ) / T1 Beta( 1 ) = Ty( 1 ) * ( Musun(Np) *Tz(1)+Rmu0 ) / T1 X( 1 ) = Alph( 1 ) + Beta( 1 ) DO 100 J = 2, Lastm2, 2 I = J / 2 Ip1 = I + 1 T1 = Rmu0**2 - Lm( Ip1 )**2 Alph( Ip1 ) = Tx( Ip1 ) / T1 Beta( Ip1 ) = Ty( Ip1 ) * ( Musun(Np) *Tz(Ip1)+Rmu0 ) / & T1 Xx = Rmu0 * Tau( Ip1 ) IF ( Xx.GT.30. ) THEN Exsun( Ip1 ) = 0. ELSE Exsun( Ip1 ) = Exp( -Rmu0*Tau(Ip1) ) END IF X( J ) = ( Alph(I)-Alph(Ip1) ) * Exsun( Ip1 ) X( J+1 ) = ( Beta(I)-Beta(Ip1) ) * Exsun( Ip1 ) 100 CONTINUE Xx = Rmu0 * Tau( Ib ) IF ( Xx.GT.30. ) THEN Exsun( Ib ) = 0. ELSE Exsun( Ib ) = Exp( -Rmu0*Tau(Ib) ) END IF X( Last ) = ( Alph(Ip1)-Beta(Ip1)+ & Alb* (Flxsun(Np)-Alph(Ip1)-Beta(Ip1)) ) * & Exsun( Ib ) C C ** Solve penta-diagonal system with r.h.s."X". C ** soln goes into "X" C CALL Leqt2b( Ss, Last, 2, 2, Nlde2, X, 1, Nlde2, Ijob, Work, & Nlde2, Work(1,8), Ier ) IF (Ier .NE. 0) GO TO 150 Dir( 1, Np ) = Flxsun( Np ) Flxd( 1, Np ) = ( 1.-Pp(1) ) / Ex( 1 ) * X( 1 ) + & ( 1.+Pp(1) ) * Ex( 1 ) * X( 2 ) - & Alph( 1 ) - Beta( 1 ) Flxu( 1, Np ) = ( 1.+Pp(1) ) / Ex( 1 ) * X( 1 ) + & ( 1.-Pp(1) ) * Ex( 1 ) * X( 2 ) - & Alph( 1 ) + Beta( 1 ) DO 110 I = 1, Nzm1 IF ( I+1.GT.Ib ) THEN Dir( I+1, Np ) = 0. Flxd( I+1, Np ) = 0. Flxu( I+1, Np ) = 0. ELSE Dir( I+1, Np ) = Flxsun( Np ) * Exsun( I+1 ) Flxd( I+1, Np ) = ( 1.-Pp(I) ) * X( 2 *I-1 ) + & ( 1.+Pp(I) ) * X( 2 *I ) - & ( Alph(I)+Beta(I) ) * Exsun( I+1 ) Flxu( I+1, Np ) = ( 1.+Pp(I) ) * X( 2 *I-1 ) + & ( 1.-Pp(I) ) * X( 2 *I ) - & ( Alph(I)-Beta(I) ) * Exsun( I+1 ) END IF 110 CONTINUE 120 CONTINUE END IF DO 130 L = 1, Nz Dirflx( L ) = Dir( L, 1 ) Difflx( L ) = Flxd( L, 1 ) Upwflx( L ) = Flxu( L, 1 ) 130 CONTINUE RETURN 150 Delta=.False. !test 120425 c print*, '-------alb' c print*, Alb c print*, '-------dttot' c print*, Dttot c print*, '-------g' c print*, Gf c print*, '-------ssa' c print*, Omega c print*, '-------solflx' c print*, Solflx c print*, '-------rmu' c print*, Rmu RETURN END C*********************************************************************** SUBROUTINE Effamt( Nzc, Zc, Pc, T, H2oden, Xh2o, Co3, Musun, & Maxlev, H2oeff, Xl, Xlstar ) C*********************************************************************** C C Calculates effective absorber amounts of water vapor and ozone. C The pressure and temperature scaling function of water vapor is C given by Chou et al., 1999: A solar radiation parameterization C (CLIRAD-SW) for atmospheric studies. NASA Tech. Memo, NASA Goddard C Space Flight Center, Greenbelt, MD, 38 pp. C INPUT: C C Co3(Lc) : Lc=1 to Nzc; C H2oden(Lc) : Lc=1 to Nzc; water vapor density at computational levels C Musun : cosine of solar zenith angle C Nzc : number of computational levels C Pc(Lc) : Lc=1 to Nzc; pressure at computational levels (mb) C T(Lc) : Lc=1 to Nzc; temperature at computational levels (K) C Zc(Lc) : Lc=1 to Nzc; altitude of computational levels (km) C C .. Scalar Arguments .. INTEGER Maxlev, Nzc REAL Musun, Xh2o C .. C .. Array Arguments .. REAL Co3( * ), H2oden( * ), H2oeff( Maxlev-1 ), Pc( * ), T( * ), & Xl( Maxlev ), Xlstar( Maxlev ), Zc( * ) C .. C .. Local Scalars .. INTEGER Lc REAL Bot, Facmag, Pt1, Pt2, Scah2o, Top, Ul, Ut C .. C .. Local Arrays .. REAL Ch2oef( Maxlev ) C .. C .. External Functions .. c REAL Averag c EXTERNAL Averag C .. C .. Intrinsic Functions .. INTRINSIC Sqrt C .. C .. Data statements .. C DATA Scah2o / 0.8 / C .. C C ** Calculate cumulative effective water vapor amount C Ch2oef( 1:Nzc ) = 0.0 DO 10 Lc = 1, Nzc - 1 Pt1 = ( Pc(Lc)/300.)**0.8 * ( 1.+0.00135*(T(Lc)-240.) ) Pt2 = ( Pc(Lc+1)/300.)**0.8 * ( 1.+0.00135*(T(Lc+1)-240.) ) Top = H2oden( Lc ) * Pt1 Bot = H2oden( Lc+1 ) * Pt2 Ch2oef( Lc+1 ) = Ch2oef( Lc ) + & ( Zc(Lc)-Zc(Lc+1) ) * 0.1 * Averag( Top, Bot ) 10 CONTINUE C ** Calculate effective water vapor amount in each layer C DO 20 Lc = 1, Nzc - 1 H2oeff( Lc ) = Xh2o * ( Ch2oef(Lc+1)-Ch2oef(Lc) ) 20 CONTINUE C ** Calculate effective ozone amount in each layer for C ** the direct and diffuse path C Facmag = 35. / ( 1224. *Musun**2.+1. )**0.5 Ut = Co3( Nzc ) DO 30 Lc = 1, Nzc Ul = Co3( Lc ) Xl( Lc ) = Ul * Facmag Xlstar( Lc ) = Ut * Facmag + 1.9 * ( Ut-Ul ) 30 CONTINUE c print*, xh2o c print*, ch2oef(nzc)*xh2o END C*********************************************************************** SUBROUTINE Leqt1b( A, N, Nlc, Nuc, Ia, B, M, Ib, Ijob, Xl, Ier ) C*********************************************************************** C C Matrix decomposition, linear equation solution C C Written by W. Wiscombe. C C C .. Scalar Arguments .. INTEGER Ia, Ib, Ier, Ijob, M, N, Nlc, Nuc C .. C .. Array Arguments .. DOUBLE PRECISION A( Ia, 1 ), B( Ib, 1 ), Xl( N, 1 ) C .. C .. Local Scalars .. INTEGER I, Ik, J, Jbeg, Jend, K, K1, Kk, L, Nc, Nlc1, Nn REAL P, Q, Rn C .. C .. Intrinsic Functions .. INTRINSIC Abs C .. Ier = 0 Jbeg = Nlc + 1 Nlc1 = Jbeg IF ( Ijob.EQ.2 ) GO TO 170 Rn = N C ** Restructure the matrix, find reciprocal of the C ** largest absolute value in row I C I = 1 Nc = Jbeg + Nuc Nn = Nc Jend = Nc IF ( N.EQ.1 .OR. Nlc.EQ.0 ) GO TO 50 10 K = 1 P = 0. DO 20 J = Jbeg, Jend A( I, K ) = A( I, J ) Q = Abs( A(I,K) ) IF ( Q.GT.P ) P = Q K = K + 1 20 CONTINUE IF ( P.EQ.0. ) GO TO 280 Xl( I, Nlc1 ) = 1. / P IF ( K.GT.Nc ) GO TO 40 DO 30 J = K, Nc A( I, J ) = 0. 30 CONTINUE 40 I = I + 1 Jbeg = Jbeg - 1 IF ( (Jend-Jbeg).EQ.N ) Jend = Jend - 1 IF ( I.LE.Nlc ) GO TO 10 Jbeg = I Nn = Jend 50 Jend = N - Nuc DO 90 I = Jbeg, N P = 0. DO 60 J = 1, Nn Q = Abs( A(I,J) ) IF ( Q.GT.P ) P = Q 60 CONTINUE IF ( P.EQ.0. ) GO TO 280 Xl( I, Nlc1 ) = 1. / P IF ( I.EQ.Jend ) GO TO 80 IF ( I.LT.Jend ) GO TO 90 K = Nn + 1 DO 70 J = K, Nc A( I, J ) = 0. 70 CONTINUE 80 Nn = Nn - 1 90 CONTINUE L = Nlc C ** L-U decomposition C DO 160 K = 1, N Cc IF(ABS(A(K,1)).LT.1.E-20) A(K,1)=0. Cc IF(ABS(XL(K,NLC1)).LT.1.E-20) XL(K,NLC1)=0. P = Abs( A(K,1) ) * Xl( K, Nlc1 ) I = K IF ( L.LT.N ) L = L + 1 K1 = K + 1 IF ( K1.GT.L ) GO TO 110 DO 100 J = K1, L Cc IF(ABS(A(J,1)).LT.1.E-20) A(J,1)=0. Cc IF(ABS(XL(J,NLC1)).LT.1.E-20) XL(J,NLC1)=0. Q = Abs( A(J,1) ) * Xl( J, Nlc1 ) IF ( Q.LE.P ) GO TO 100 P = Q I = J 100 CONTINUE 110 Xl( I, Nlc1 ) = Xl( K, Nlc1 ) Xl( K, Nlc1 ) = I C C ** Singularity found C IF ( (Rn+P).EQ.Rn ) GO TO 280 C C ** Interchange rows I and K C IF ( K.EQ.I ) GO TO 130 DO 120 J = 1, Nc P = A( K, J ) A( K, J ) = A( I, J ) A( I, J ) = P 120 CONTINUE 130 IF ( K1.GT.L ) GO TO 160 DO 150 I = K1, L P = A( I, 1 ) / A( K, 1 ) Ik = I - K Xl( K1, Ik ) = P DO 140 J = 2, Nc Cc IF(ABS(P).LT.1.E-20) P=0. Cc IF(ABS(A(K,J)).LT.1.E-20) A(K,J)=0. A( I, J-1 ) = A( I, J ) - P * A( K, J ) 140 CONTINUE A( I, Nc ) = 0. 150 CONTINUE 160 CONTINUE IF ( Ijob.EQ.1 ) GO TO 270 170 L = Nlc DO 220 K = 1, N I = Xl( K, Nlc1 ) IF ( I.EQ.K ) GO TO 190 DO 180 J = 1, M P = B( K, J ) B( K, J ) = B( I, J ) B( I, J ) = P 180 CONTINUE 190 IF ( L.LT.N ) L = L + 1 K1 = K + 1 IF ( K1.GT.L ) GO TO 220 DO 210 I = K1, L Ik = I - K P = Xl( K1, Ik ) DO 200 J = 1, M Cc IF(ABS(P).LT.1.E-20) P=0. Cc IF(ABS(B(K,J)).LT.1.E-20) B(K,J)=0. B( I, J ) = B( I, J ) - P * B( K, J ) 200 CONTINUE 210 CONTINUE 220 CONTINUE Jbeg = Nuc + Nlc DO 260 J = 1, M L = 1 K1 = N + 1 DO 250 I = 1, N K = K1 - I P = B( K, J ) IF ( L.EQ.1 ) GO TO 240 DO 230 Kk = 2, L Ik = Kk + K Cc IF(ABS(A(K,KK)).LT.1.E-20) A(K,KK)=0. Cc IF(ABS(B(IK-1,J)).LT.1.E-20) B(IK-1,J)=0. P = P - A( K, Kk ) * B( Ik-1, J ) 230 CONTINUE 240 B( K, J ) = P / A( K, 1 ) IF ( L.LE.Jbeg ) L = L + 1 250 CONTINUE 260 CONTINUE 270 RETURN 280 Ier = 129 WRITE ( *, FMT=9000 ) return 9000 FORMAT ( /, /, ' LEQT1B ERROR--BANDED MATRIX IS SINGULAR' ) END C*********************************************************************** SUBROUTINE Leqt2b( A, N, Nlc, Nuc, Ia, B, M, Ib, Ijob, U, Iu, Xl, & Ier ) C*********************************************************************** C Written by W. Wiscombe C .. Scalar Arguments .. INTEGER Ia, Ib, Ier, Ijob, Iu, M, N, Nlc, Nuc C .. C .. Array Arguments .. DOUBLE PRECISION A( Ia, 1 ), B( Ib, 1 ), U( Iu, 1 ), Xl( N, 1 ) C .. C .. Local Scalars .. INTEGER I, Ir, Iter, Itmax, J, Jj, Jk, K, Kk, L, Nc, Ncc, Nlcp1, & Nlcp2, Nmnuc, Nu2, Nu3 REAL Dxnorm, Prec, Xnorm, Xxxx DOUBLE PRECISION Sum C .. C .. External Subroutines .. c EXTERNAL Leqt1b C .. C .. Intrinsic Functions .. INTRINSIC Abs, Epsilon, Max C .. C .. Data statements .. DATA Itmax / 200 / C .. Prec = Epsilon( 1.0 ) Ier = 0 Nlcp1 = Nlc + 1 Ncc = Nuc + Nlcp1 Nlcp2 = Nlc + 2 Nmnuc = N - Nuc Nu2 = Ncc + 1 Nu3 = Ncc + 2 IF ( Ijob.EQ.2 ) GO TO 30 DO 20 I = 1, N DO 10 J = 1, Ncc U( I, J ) = A( I, J ) 10 CONTINUE 20 CONTINUE CALL Leqt1b( U, N, Nlc, Nuc, Iu, B, M, Ib, 1, Xl, Ier ) IF (ier .ne. 0 ) GO TO 140 IF ( Ijob.EQ.1 ) GO TO 130 30 DO 120 J = 1, M DO 40 I = 1, N U( I, Nu2 ) = B( I, J ) 40 CONTINUE CALL Leqt1b( U, N, Nlc, Nuc, Iu, U(1,Nu2), 1, Iu, 2, Xl, Ier ) IF (ier .ne. 0 ) GO TO 140 Xnorm = 0. DO 50 I = 1, N IF ( Xnorm.LT.Abs(U(I,Nu2)) ) THEN Xnorm = Abs( U(I,Nu2) ) END IF 50 CONTINUE IF ( Xnorm.EQ.0. ) GO TO 120 DO 90 Iter = 1, Itmax Nc = Ncc Kk = 1 DO 70 I = 1, N Sum = B( I, J ) L = Nlcp2 - I Ir = Max( L, 1 ) IF ( L.LE.0 ) Kk = Kk + 1 K = Kk DO 60 Jj = Ir, Nc Sum = Sum - A( I, Jj ) * U( K, Nu2 ) K = K + 1 60 CONTINUE U( I, Nu3 ) = Sum IF ( I.GE.Nmnuc ) Nc = Nc - 1 70 CONTINUE CALL Leqt1b( U, N, Nlc, Nuc, Iu, U(1,Nu3), 1, Iu, 2, Xl, & Ier ) IF (ier .ne. 0 ) GO TO 140 Dxnorm = 0. DO 80 I = 1, N U( I, Nu2 ) = U( I, Nu2 ) + U( I, Nu3 ) IF ( Dxnorm.LT.Abs(U(I,Nu3)) ) THEN Dxnorm = Abs( U(I,Nu3) ) END IF 80 CONTINUE Xxxx = Xnorm + Dxnorm IF ( Abs(Xxxx-Xnorm).LE.Prec ) GO TO 100 90 CONTINUE Ier = 130 100 DO 110 Jk = 1, N B( Jk, J ) = U( Jk, Nu2 ) 110 CONTINUE IF ( Ier.NE.0 ) GO TO 140 120 CONTINUE 130 RETURN 140 WRITE ( *, FMT=9000 ) Return 9000 FORMAT ( /, /, & 'LEQT2B ERROR--ITERATIVE IMPROVEMENT FAILED TO CONVER', & 'GE BECAUSE BANDED MATRIX IS TOO ILL-CONDITIONED' ) END C*********************************************************************** SUBROUTINE Main2( Nab, Nkn, Pkn, Kn, H2oeff, Nzc, Dtauex, Dtausc, & Gf, Maxlev, Musun, Alb, Solflx, Flxsun, Fld, & Fldn, Flup, Delta) C*********************************************************************** C C Completes the calculation of fluxes for current spectral C interval, including forming the weighted sums of C monochromatic multiple-scattering calculations according to C the k-distribution scheme. C C Dta(Lc) : Lc=1 to Nzc-1; pseudo-molecular optical depth C in layer -Lc- C Dtausc(Lc): Lc=1 to Nzc-1; continuum scattering optical depth C for layer -Lc- C Dtauex(Lc): Lc=1 to Nzc-1; extiction optical depth for layer -Lc- C Fld(Lc) : Lc=1 to Nzc; downward direct flux at level -Lc-, for C current spectral interval (watts/m**2) C Fldn(Lc) : Lc=1 to Nzc; downward flux (diffuse + direct) at level C -Lc-, for current spectral interval (watts/m**2) C Flup(Lc) : Lc=1 to Nzc; same as *Fldn*, only upward C Flxsun : solar flux normal to the atmosphere (watts/m**2) C Nab : number of active absorber classes in current spectral C interval (at most 1 -- H2O) C Tau(Lc) : Lc=1 to Nzc; extinction optical depth between the top C of the atmosphere and level -Lc- C C .. Scalar Arguments .. INTEGER Maxlev, Nab, Nkn, Nzc REAL Alb, Flxsun, Musun, Solflx C .. C .. Array Arguments .. REAL Dtauex( * ), Dtausc( * ), Fld( Maxlev ), Fldn( Maxlev ), & Flup( Maxlev ), Gf( * ), H2oeff( Maxlev-1 ), Kn( * ), & Pkn( * ) C .. C .. Local Scalars .. INTEGER Lc, N C .. C .. Local Arrays .. REAL Dir( Maxlev ), Dta( Maxlev-1 ), Flxd( Maxlev ), & Flxu( Maxlev ), Tau( Maxlev ) DOUBLE PRECISION Asymf( Maxlev-1 ), Dttot( Maxlev-1 ), & Omega( Maxlev-1 ) LOGICAL delta C .. C .. External Subroutines .. c EXTERNAL Delted C .. C .. Intrinsic Functions .. INTRINSIC Exp C .. !test 120425 c DOUBLE PRECISION Oo( Maxlev-1 ) !test 120425 end Dttot = 0.0 Omega = 0.0 Asymf = 0.0 Tau = 0.0 IF ( Nab.EQ.0 ) THEN C ** No molecular absorption--solve only a single problem C DO 10 Lc = 1, Nzc - 1 Dttot( Lc ) = Dtauex( Lc ) Omega( Lc ) = Dtausc( Lc ) / Dttot( Lc ) Asymf( Lc ) = Gf( Lc ) Tau( Lc+1 ) = Tau( Lc ) + Dttot( Lc ) 10 CONTINUE CALL Delted( Solflx, Dttot, Omega, Asymf, Maxlev, Alb, Musun, & Nzc, Fld, Fldn, Flup, Delta ) IF ( .NOT. Delta) GO TO 100 C ** Add direct beam to diffuse down-flux C DO 20 Lc = 1, Nzc Fldn( Lc ) = Fldn( Lc ) + Fld( Lc ) 20 CONTINUE C ** Calculate direct beam from Beer-Lambert law C DO 30 Lc = 1, Nzc Fld( Lc ) = Flxsun * Exp( -Tau(Lc)/Musun ) 30 CONTINUE ELSE C ** Molecular absorption--solve multiple problems C Fld = 0.0 Fldn = 0.0 Flup = 0.0 DO 90 N = 1, Nkn C C ** Pseudo-absorption-optical-depth for active absorber C DO 40 Lc = 1, Nzc - 1 Dta( Lc ) = Kn( N ) * H2oeff( Lc ) 40 CONTINUE DO 50 Lc = 1, Nzc - 1 Dttot( Lc ) = Dtauex( Lc ) + Dta( Lc ) Omega( Lc ) = Dtausc( Lc ) / Dttot( Lc ) Asymf( Lc ) = Gf( Lc ) Tau( Lc+1 ) = Tau( Lc ) + Dttot( Lc ) 50 CONTINUE CALL Delted( Solflx, Dttot, Omega, Asymf, Maxlev, Alb, & Musun, Nzc, Dir, Flxd, Flxu, delta ) IF ( .NOT. Delta) GO TO 100 C ** Add direct beam to diffuse down-flux C DO 60 Lc = 1, Nzc Flxd( Lc ) = Flxd( Lc ) + Dir( Lc ) 60 CONTINUE C ** Calculate direct beam from Beer-Lambert law C c write(*,*)tau(nzc),musun DO 70 Lc = 1, Nzc Dir( Lc ) = Flxsun * Exp( -Tau(Lc)/Musun ) 70 CONTINUE DO 80 Lc = 1, Nzc Fld( Lc ) = Fld( Lc ) + Pkn( N ) * Dir( Lc ) Fldn( Lc ) = Fldn( Lc ) + Pkn( N ) * Flxd( Lc ) Flup( Lc ) = Flup( Lc ) + Pkn( N ) * Flxu( Lc ) 80 CONTINUE 90 CONTINUE END IF 100 RETURN END C*********************************************************************** SUBROUTINE Prnflx( Flxdir, Flxdnw, Flxupw, Musun, Nz, P, Maxlev, & Wllo, Wlhi ) C*********************************************************************** C C Prints vertical fluxes and quantities derived from them C C .. Scalar Arguments .. INTEGER Maxlev, Nz REAL Musun, Wlhi, Wllo C .. C .. Array Arguments .. REAL Flxdir( * ), Flxdnw( * ), Flxupw( * ), P( * ) C .. C .. Local Scalars .. CHARACTER Lab * 21 INTEGER L REAL Hrat C .. C .. Local Arrays .. REAL Absycm( Maxlev ), Absyly( Maxlev ), Albedo( Maxlev ), & Fnet( Maxlev ), Heatng( Maxlev ), Trancm( Maxlev ), & Tranly( Maxlev ) C .. C ** Create print header C WRITE ( Lab, FMT=9000 ) Wllo, Wlhi WRITE ( *, FMT=9010 ) Lab, Musun C ** Calculate net fluxes, albedos, transmissivities and C ** absorptivities at all computational levels C Albedo = 0.0 DO 10 L = 1, Nz Fnet( L ) = Flxdnw( L ) - Flxupw( L ) IF ( Flxdnw(L).GT.0. ) Albedo( L ) = Flxupw( L ) / Flxdnw( L ) IF ( Flxdnw(1).EQ.0. ) THEN Trancm( L ) = 0.0 Absycm( L ) = 0.0 ELSE Trancm( L ) = Flxdnw( L ) / Flxdnw( 1 ) Absycm( L ) = ( Fnet(1)-Fnet(L) ) / Flxdnw( 1 ) END IF 10 CONTINUE C C ** Calculate heating rate, transmissivity and C ** absorptivity of layers C DO 20 L = 1, Nz - 1 Heatng( L ) = -8.39 * ( Fnet(L+1)-Fnet(L) ) / ( P(L+1)-P(L) ) IF ( Flxdnw(L).EQ.0. ) THEN Tranly( L ) = 0.0 Absyly( L ) = 0.0 ELSE Tranly( L ) = Flxdnw( L+1 ) / Flxdnw( L ) Absyly( L ) = ( Fnet(L)-Fnet(L+1) ) / Flxdnw( L ) END IF WRITE ( *, FMT=9020 ) L, Flxdir( L ), & Flxdnw( L ) - Flxdir( L ), Flxdnw( L ), Flxupw( L ), & Fnet( L ), Albedo( L ), Trancm( L ), Absycm( L ), & L, L + 1, Heatng( L ), Tranly( L ), Absyly( L ) 20 CONTINUE Hrat = -8.39 * ( Fnet(Nz)-Fnet(1) ) / ( P(Nz)-P(1) ) WRITE ( *, FMT=9020 ) Nz, Flxdir( Nz ), & Flxdnw( Nz ) - Flxdir( Nz ), Flxdnw( Nz ), Flxupw( Nz ), & Fnet( Nz ), Albedo( Nz ), Trancm( Nz ), Absycm( Nz ) WRITE ( *, FMT=9030 ) Hrat 9000 FORMAT ( F8.2, '-', F8.2, ' mcr' ) 9010 FORMAT ( /, /, /, /, ' <', 44 ('-'), 5X, A, 5X, 45 ('-'), '>', /, & /, ' COS OF SOL ZEN ANGLE =', F7.4, /, & ' LEVEL DOWNWARD DOWNWARD DOWNWARD UPWARD NET FLUX' & , ' ALBEDO CUMUL CUMUL LAYERS HEATING', & ' LAYER LAYER', /, & ' DIRECT DIFFUSE TOTAL DIFFUSE ' & , ' TRANSY ABSORY RATE', & ' TRANSY ABSORY' ) 9020 FORMAT ( I6, 5F11.4, 3F8.4, I4, ' -', I3, 1P, E13.4, 0P, 2F8.4 ) 9030 FORMAT ( /, ' TOTAL COLUMN HEATING RATE (DEG/DAY):', 1P, E13.4 ) END C*********************************************************************** REAL FUNCTION Scaray( Iwl, Qvtop, Qvbot, Ztop, Zbot, Ptop, & Pbot ) C*********************************************************************** C C Computes Rayleigh scattering optical depth (including C corrections for humidity and the variation of gravitation C with altitude) C .. Parameters .. INTEGER Maxwvl PARAMETER ( Maxwvl=7 ) C .. C .. Scalar Arguments .. INTEGER Iwl REAL Pbot, Ptop, Qvbot, Qvtop, Zbot, Ztop C .. C .. Local Scalars .. REAL Humidy C .. C .. Local Arrays .. REAL Tmp( Maxwvl ) C .. C .. Data statements .. DATA Tmp / 2.689E-5, 8.220E-6, 3.640E-6, 1.846E-6, 3.630E-7, & 3.630E-7,3.630E-7 / C .. Humidy = 0.5 * ( Qvtop+Qvbot ) Scaray = 27.019 * Tmp( Iwl ) * ( 1.0+0.5 * (Ztop+Zbot)/6371. )** & 2 * ( 1.+29./18. *Humidy ) / ( 1.+Humidy ) * & ( Pbot-Ptop ) END C*********************************************************************** SUBROUTINE Scat( Iwl, Nzc, Zc, Pc, Qv, Noray, Noaer, & Nocld, Maxlev, Dteaer, Omaer, Gfaer, Dtecld, & Omcld, Gfcld, Dtauex, Dtausc, Gf ) C*********************************************************************** C C For current spectral interval, calculates extinction and C scattering optical depths, and asymmetry factors by C combining Rayleigh and particulate (aerosol and cloud) C extinction C C .. Scalar Arguments .. LOGICAL Noaer, Nocld, Noray INTEGER Iwl, Maxlev, Nzc C .. C .. Array Arguments .. REAL Dteaer( * ), Dtecld( * ), Gfaer( * ), Gfcld( * ), & Omaer( * ), Omcld( * ), Pc( * ), Qv( * ), Zc( * ), & Dtauex( Maxlev-1 ), Dtausc( Maxlev-1 ), Gf( Maxlev-1 ) C .. C .. Local Scalars .. INTEGER Lc REAL Dtsaer, Dtscld, Ssalb, Totray C .. Dtauex = 0.0 Dtausc = 0.0 Gf = 0.0 C ** Avoid possible divide-by-zero C Dtauex = Max( Dtauex, Tiny(Dtauex) ) IF ( .NOT.Noray ) THEN C C ** Compute Rayleigh scattering optical depth C Totray = 0.0 DO 10 Lc = 1, Nzc - 1 Dtausc( Lc ) = Scaray( Iwl, Qv(Lc), Qv(Lc+1), Zc(Lc), & Zc(Lc+1), Pc(Lc), Pc(Lc+1) ) Dtauex( Lc ) = Dtausc( Lc ) Dtauex( Lc ) = Max( Dtauex(Lc), Tiny(Dtauex(Lc)) ) Totray = Totray + Dtauex( Lc ) 10 CONTINUE !test c print*, '--- ray' c print*, Dtausc( 1:Nzc-1 ) !test end END IF IF ( .NOT.Noaer ) THEN C ** Add contribution of particle extinction C DO 20 Lc = 1, Nzc - 1 Dtsaer = Omaer( Lc ) * Dteaer( Lc ) IF ( (Dtausc(Lc)+Dtsaer).GT.0.0 ) & Gf( Lc ) = ( Dtausc(Lc)*Gf(Lc) + Dtsaer*Gfaer(Lc) ) / & ( Dtausc(Lc)+Dtsaer ) Dtausc( Lc ) = Dtausc( Lc ) + Dtsaer Dtauex( Lc ) = Dtauex( Lc ) + Dteaer( Lc ) Dtauex( Lc ) = Max( Dtauex(Lc), Tiny(Dtauex(Lc)) ) cc Ssalb = Dtausc( Lc ) / Dtauex( Lc ) cc print *, Lc, Dtauex(Lc), Ssalb, Gf(Lc) cc print *, Lc, Dtauex(Lc), Dtausc(Lc), Gf(Lc) !test c print*, Lc, Dtsaer, Omaer( Lc ), Dteaer( Lc ) 20 CONTINUE !test c print*, '--- aer' c print*, Dtausc( 1:Nzc-1 ) c print*, '--- aer' c print*, Dtausc( 1:Nzc-1 ) !test end !test 120425 c if (any(dtauex<0) .or. any(dtausc<0)) then c print*, '------omaer' c print*, omaer(1:Nzc-1) c print*, '------dteaer' c print*, dteaer(1:nzc-1) c print*, '------gf' c print*, gf(1:nzc-1) c endif !test 120425 end END IF IF ( .NOT.Nocld ) THEN C ** Add contribution of cloud extinction C DO 30 Lc = 1, Nzc - 1 Dtscld = Omcld( Lc ) * Dtecld( Lc ) IF ( (Dtausc(Lc)+Dtscld).GT.0.0 ) & Gf( Lc ) = ( Dtausc(Lc)*Gf(Lc) + Dtscld*Gfcld(Lc) ) / & ( Dtausc(Lc)+Dtscld ) Dtausc( Lc ) = Dtausc( Lc ) + Dtscld Dtauex( Lc ) = Dtauex( Lc ) + Dtecld( Lc ) Dtauex( Lc ) = Max( Dtauex(Lc), Tiny(Dtauex(Lc)) ) c Ssalb = Dtausc( Lc ) / Dtauex( Lc ) c print *, Lc, Dtauex(Lc), Ssalb, Gf(Lc) c print *, Lc, Dtauex(Lc), Dtausc(Lc), Gf(Lc) !test c print*, '--- cld' c print*, Dtausc( 1:Nzc-1 ) !test end 30 CONTINUE !test 120425 c if (any(dtauex<0) .or. any(dtausc<0)) then c print*, '------omcld' c print*, omcld(1:Nzc-1) c print*, '------dtecld' c print*, dtecld(1:nzc-1) c print*, '------gf' c print*, gf(1:nzc-1) c endif !test 120425 end END IF END !----------------------------------------------------------------------- C C Calculates extinction optical depth, single scattering C albedo, and asymmetry factor for particulate C scattering by interpolating layer values of aerosol models C C Tauaer : total column aerosol optical depth at 0.55 microns C Laer(Lc): Lc=1 to Nzc, model-aerosol layer to be used in C computational layer Lc C Mlsra : max number of layers in standard radiation atmospheres C Mzsra : max number of levels in standard radiation atmospheres C Mnwl : max number of spectral intervals in aerosol models C Msra : max number of aerosol models !----------------------------------------------------------------------- SUBROUTINE Setaer( Nzc, Zc, Iaer, aerObj, Laer, Nwl, Maxlev, & Maxwvl, Dtauex, Ssalb, Gf ) !yma 160410 !----------------------------------------------------------------------- USE Standard_Aerosol ,ONLY: Bsra, Osra, Gsra, Sclht !yma 160410 IMPLICIT NONE integer ,intent(in) :: Nzc real ,intent(in) :: Zc( * ) integer ,intent(in) :: Iaer class(*) ,intent(in) :: aerObj integer ,intent(in) :: Laer( Maxlev-1 ) integer ,intent(in) :: Nwl integer ,intent(in) :: Maxlev integer ,intent(in) :: Maxwvl real ,intent(out) :: Dtauex( Maxlev-1, Maxwvl ) real ,intent(out) :: Ssalb( Maxlev-1, Maxwvl ) real ,intent(out) :: Gf( Maxlev-1, Maxwvl ) C .. Local variables INTEGER :: Iwl, Lc REAL :: Xtau,Stau real :: ext, ssa, asy, Hscl !yma 160410 real :: userExtTau !total optical depth real :: userAbsTau !user input total absorption optical depth REAL :: Tauext( Maxlev, Maxwvl ) real :: TauAbs( Maxlev, Maxwvl ) real :: dTauAb( Maxlev-1, Maxwvl ) Tauext(:,:) = 0.0 TauAbs(:,:) = 0.0 Dtauex(:,:) = 0.0 Ssalb(:,:) = 0.0 Gf(:,:) = 0.0 DO 20 Lc = 1, Nzc - 1 IF ( Laer(Lc).NE.0 ) THEN DO 10 Iwl = 1, Nwl C ** Integrate aerosol extinction and scattering coeffs C ** between adjacent levels to get layer optical depths, C ** and choose aerosol model layer single scattering albedos C ** and asymmetry factors to get computation layer values c------------------------->>> ! 160410 select type(aerObj) type is (real) !use standard aer model ext = Bsra( Iwl, Laer(Lc), Iaer ) ssa = Osra( Iwl, Laer(Lc), Iaer ) asy = Gsra( Iwl, Laer(Lc), Iaer ) Hscl = Sclht( Laer(Lc), Iaer ) type is (SRB_UserAer_Form1) !user provided aer model in form 1 ext = aerObj%Bsra( Iwl, Laer(Lc) ) ssa = aerObj%Osra( Iwl, Laer(Lc) ) asy = aerObj%Gsra( Iwl, Laer(Lc) ) Hscl = aerObj%Sclht( Laer(Lc) ) type is (SRB_UserAer_Form2) !use standard aer model ext = Bsra( Iwl, Laer(Lc), Iaer ) ssa = Osra( Iwl, Laer(Lc), Iaer ) asy = Gsra( Iwl, Laer(Lc), Iaer ) Hscl = Sclht( Laer(Lc), Iaer ) end select Ssalb( Lc, Iwl ) = ssa Gf( Lc, Iwl ) = asy Dtauex( Lc, Iwl ) = Hscl * ext * & ( Exp(-Zc(Lc+1)/Hscl)- Exp(-Zc(Lc)/Hscl)) dTauAb( Lc, Iwl ) = Dtauex( Lc, Iwl ) * (1-ssa) Tauext( Lc+1, Iwl ) = Tauext( Lc, Iwl ) + & Dtauex( Lc, Iwl ) TauAbs( Lc+1, Iwl ) = TauAbs( LC, Iwl ) + & dTauAb( Lc, Iwl ) c------------------------->>><<< ! 160410 !test c print *, Lc, Iwl, Laer(Lc), Dtauex(Lc,Iwl), Ssalb(Lc,Iwl), c & Gf(Lc,Iwl) c print *, Lc, Iwl, Laer(Lc), Dtauex(Lc,Iwl), Hscl, ext, c & ( Exp(-Zc(Lc+1)/Hscl)- Exp(-Zc(Lc)/Hscl)) 10 CONTINUE END IF 20 CONTINUE !---------------------->>> userAer select type(aerObj) type is (real) userExtTau = aerObj userAbsTau = 0. type is (SRB_UserAer_Form1) userExtTau = aerObj%Tauaer userAbsTau = 0. type is (SRB_UserAer_Form2) userExtTau = aerObj%Tauaer userAbsTau = aerObj%absTau end select !---------------------->>><<< userAer IF ( userExtTau.GE.0.0 ) THEN DO Lc = 1, Nzc - 1 IF ( Laer(Lc).NE.0 ) THEN Xtau = userExtTau / Tauext( Nzc, 3 ) DO Iwl = 1, Nwl Dtauex( Lc, Iwl ) = Xtau * Dtauex( Lc, Iwl ) Tauext( Lc+1, Iwl ) = Xtau * Tauext( Lc+1, Iwl ) ENDDO ENDIF !test c print*, Lc,Iwl,Dtauex(Lc,Iwl),Tauext( Lc+1, Iwl ),Xtau ENDDO END IF !---------------------->>> userAer IF ( userAbsTau.GE.0.0 ) THEN DO Lc = 1, Nzc - 1 IF ( Laer(Lc).NE.0 ) THEN Stau = userAbsTau / TauAbs( Nzc, 3 ) DO Iwl = 1, Nwl Ssalb( Lc, Iwl ) = 1.-Stau*dTauAb( Lc,Iwl )/ & Dtauex( Lc,Iwl) TauAbs( Lc+1, Iwl ) = Stau*TauAbs( Lc+1, Iwl ) ENDDO ENDIF !test c print*, Lc,Iwl,Ssalb(Lc,Iwl),Stau,dTauAb(Lc,Iwl),Dtauex(Lc,Iwl) ENDDO END IF !---------------------->>><<< userAer c DO Lc = 1, Nzc-1 c print *, Lc, Dtauex(Lc,1), Ssalb(Lc,1), Gf(Lc,1), Laer(Lc) c print *, Lc, Dtauex(Lc,2), Ssalb(Lc,2), Gf(Lc,2), Laer(Lc) c print *, Lc, Dtauex(Lc,3), Ssalb(Lc,3), Gf(Lc,3), Laer(Lc) c END DO c print *, Nzc, Tauext(Nzc,3) c END C*********************************************************************** SUBROUTINE Setatm( Prnt, Iatm, Iaer, Maxlev, Ozone, & Pwater, Zsrf, Psrf, Nout, Pout, & Nlcld, Zcldt, Zcldb, Nzc, Zc, Pc, T, Qv, & H2oden, Co3, Xh2o, Laer, Lcld, Lprn, aerObj ) !yma 160410 C*********************************************************************** C C Sets up altitude, pressure, temperature, H2O and O3 vertical C profiles C C NOTE: Precip water values (Dh2o) are assumed to be from top to C bottom. ISCCP gives them from bottom to top, so they must be C reversed before calling this routine. C Nout : number of user-defined output pressure levels C Pout(Lu) : Lu=1 to Nout; user pressures (mb); fluxes are to be C returned at these levels; note: actual output levels C at the top and at the surface may differ C Lprn(Lu) : Lu=1 to Nzprn; actual computational level numbers at C which vertical printing is to be done C Pwater : user-defined total column amount of H2O ( g/sq cm) C Zcldb(L) : L=1 to Nlcld; cloud-base height for cloud layer -L- (km) C Zcldt(L) : L=1 to Nlcld; cloud-top height for cloud layer -L- (km) C C C .. C .. Use Statements .. USE McClatchey_Atmospheres USE Standard_Aerosol ,ONLY: Nzsra, Za ! 160410 IMPLICIT None class(*) ,intent(in) :: aerObj ! 160410 C .. C .. Scalar Arguments .. LOGICAL Prnt INTEGER Iaer, Iatm, Maxlev, Nlcld, Nout, Nzc REAL Ozone, Pwater, Psrf, Zsrf C .. C .. Array Arguments .. INTEGER Laer( Maxlev-1 ), Lcld( Maxlev-1 ), Lprn( Maxlev ) REAL Co3( * ), H2oden( * ), Pc( * ), & Pout( * ), Qv( * ), T( * ), Zc( * ), & Zcldt( * ), Zcldb( * ) C .. C .. Local Scalars .. LOGICAL Dimerr INTEGER K, L, L1p, L2p, Lbot, Lc, Ltop, Lu, Nlaer, Nz, Nzprn & REAL Frac, H2ostd, Misval, Pscale, Xo3, Xh2o C .. C .. Local Arrays .. REAL Ch2o( Maxlev ), H2oamt( Maxlev-1 ), O3amt( Maxlev-1 ), & O3den( Maxlev ), P( Maxlev ), & Pprn( Maxlev ),Z( Maxlev ),Zaerb( Maxlev ),Zaert( Maxlev), !Zaerb( Mlsra ), Zaert( Mlsra ), ! 160410 & Zprn( Maxlev ) C .. C .. External Functions .. c REAL Averag, Densty c EXTERNAL Averag, Densty C .. C .. Intrinsic Functions .. INTRINSIC Log10, Max, Min C Misval = -9999.0 P = 0.0 Z = 0.0 Pc( 1:Maxlev ) = 0.0 Zc( 1:Maxlev ) = 0.0 Nz = 0 Nzc = 0 Pprn = 0.0 C ** Set computational heights and pressures to those of C ** the standard (McClatchey) atmospheres above height C ** -Zsrf- Zsrf = MAX( Zsrf, Zs(Mccl) ) L = 1 DO WHILE ( Zs(L).GT.Zsrf .AND. L.LT.Mccl ) Z( L ) = Zs( L ) P( L ) = Ps( L, Iatm ) L = L + 1 END DO Nz = L Dimerr = .FALSE. IF ( Nz.GT.Maxlev ) Dimerr = Wrtdim( 'Maxlev', Nz ) IF ( Dimerr ) CALL Errmsg( 'Setatm--Dimension error', .TRUE. ) C ** Set bottom computational height to height -Zsrf-, C ** and obtain pressure at that height by interpolating C ** pressure values of the standard (McClatchey) C ** atmospheres Z( Nz ) = Zsrf Frac = ( Z(Nz)-Zs(Nz-1) ) / ( Zs(Nz)-Zs(Nz-1) ) P( Nz ) = Ps( Nz-1, Iatm ) * ( Ps(Nz,Iatm)/Ps(Nz-1,Iatm) )**Frac C C ** Scale pressure values at all levels so that the C ** pressure at bottom is the same as Psrf IF (Psrf .EQ. Misval) THEN Pscale=1.0 ELSE Pscale = Psrf / P( Nz ) ENDIF DO L = 1, Nz P( L ) = Pscale * P( L ) END DO Nzprn = 0 IF ( Nout.GT.0 ) THEN C ** Find user-defined output levels included in C ** current pressure profile IF ( Pout(Nout).LT.P(1) .OR. Pout(1).GT.P(Nz) ) THEN CALL Errmsg( 'User levels are out of profile', .FALSE. ) Nzprn = 0 ELSE C ** Find first output pressure level that is larger than C ** the pressure of the first computational level C L = 1 DO WHILE ( Pout(L).LE.P(1) ) L = L + 1 END DO Ltop = L C ** Find last output pressure level that is smaller than C ** the pressure of the last computational level C L = Nout DO WHILE ( Pout(L).GE.P(Nz) ) L = L - 1 END DO Lbot = L C ** Set up print levels at top of atmosphere, at output C ** pressures and at surface C L = 1 Pprn( L ) = Min( P(1), Pout(Ltop) ) DO K = Ltop, Lbot L = L + 1 Pprn( L ) = Pout(K ) END DO Nzprn = L + 1 Pprn( Nzprn ) = Max( P(Nz), Pout(Lbot) ) C ** Find altitudes for print levels C CALL Alt( Nzprn, Pprn, Z, P, Nz, Maxlev, Zprn ) END IF END IF C ** Set top and bottom altitudes of standard aerosol C ** layers C !---------------------->>> ! 160410 select type(aerObj) type is (real) !standard built-in aer model Nlaer = Nzsra( Iaer ) - 1 DO L = 1, Nlaer Zaert( L ) = Za( L, Iaer ) Zaerb( L ) = Za( L+1, Iaer ) END DO type is (SRB_UserAer_Form1) !user provided aer model form 1 Nlaer = aerObj%Nzsra - 1 DO L = 1, Nlaer Zaert( L ) = aerObj%Za( L ) Zaerb( L ) = aerObj%Za( L+1 ) END DO type is (SRB_UserAer_Form2) !user aer form 2 Nlaer = Nzsra( Iaer ) - 1 DO L = 1, Nlaer Zaert( L ) = Za( L, Iaer ) Zaerb( L ) = Za( L+1, Iaer ) END DO end select !---------------------->>><<< ! 160410 C ** Insert boundaries of aerosol, cloud and water C ** vapor layers, and output levels C Nzc = Nz Zc( 1:Maxlev ) = 0.0 Zc( 1:Nzc ) = Z( 1:Nz ) IF ( Nlaer.GT.0 ) CALL Addlev( Nzc, Zc, Nlaer, Zaert, Maxlev ) IF ( Nlaer.GT.0 ) CALL Addlev( Nzc, Zc, Nlaer, Zaerb, Maxlev ) IF ( Nlcld.GT.0 ) CALL Addlev( Nzc, Zc, Nlcld, Zcldt, Maxlev ) IF ( Nlcld.GT.0 ) CALL Addlev( Nzc, Zc, Nlcld, Zcldb, Maxlev ) IF ( Nzprn.GT.0 ) CALL Addlev( Nzc, Zc, Nzprn, Zprn, Maxlev ) C ** Calculate pressure at altitude Zc by interpolating C ** pressures P at heights Z C DO Lc = 1, Nzc L = 2 DO WHILE ( Z(L).GT.Zc(Lc) .AND. L.LT.Nz ) L = L + 1 END DO Frac = ( Zc(Lc)-Z(L-1) ) / ( Z(L)-Z(L-1) ) Pc( Lc ) = P( L-1 ) * ( P(L)/P(L-1) )**Frac END DO C ** Find levels at which to print C Lprn = 0 DO Lu = 1, Nzprn DO Lc = 1, Nzc IF ( Zc(Lc).EQ.Zprn(Lu) ) THEN Lprn( Lu ) = Lc EXIT END IF END DO cc print *, Lu, Pprn(Lu), Zprn(Lu), Lprn(Lu) END DO C ** Obtain temperature, ozone and water vapor density, C ** and relative humidity at levels -Zc- by interpolating C ** values of the standard (McClatchey) atmospheres DO Lc = 1, Nzc L = 2 DO WHILE ( Zs(L).GT.Zc(Lc) .AND. L.LT.Mccl ) L = L + 1 END DO Frac = ( Zc(Lc)-Zs(L-1) ) / ( Zs(L)-Zs(L-1) ) T( Lc ) = Ts( L-1,Iatm ) + Frac * ( Ts(L,Iatm)-Ts(L-1,Iatm) ) O3den( Lc ) = 46.6667 * O3( L-1,Iatm ) * & ( O3(L,Iatm)/O3(L-1,Iatm) )**Frac IF ( H2o(L-1,Iatm).GE.1.E-10 .AND. H2o(L,Iatm).GE.1.E-10 ) THEN H2oden( Lc ) = H2o( L-1,Iatm ) * & ( H2o(L,Iatm)/H2o(L-1,Iatm) )**Frac ELSE H2oden( Lc ) = ( H2o(L-1,Iatm)+ & Frac * (H2o(L,Iatm)-H2o(L-1,Iatm)) ) END IF Qv( Lc ) = 1.E-3 * H2oden( Lc ) / Densty( Pc(Lc), T(Lc), 0. ) END DO C ** Calculate cumulative water vapor and ozone amount of C ** standard atmosphere C Ch2o( 1:Nzc ) = 0.0 Co3( 1:Nzc ) = 0.0 DO Lc = 1, Nzc-1 Ch2o( Lc+1) = Ch2o( Lc ) + ( Zc(Lc)-Zc(Lc+1) ) * 0.1 * & Averag( H2oden(Lc), H2oden(Lc+1) ) Co3( Lc+1 ) = Co3( Lc ) + ( Zc(Lc)-Zc(Lc+1) ) * & Averag( O3den(Lc), O3den(Lc+1) ) END DO C ** Obtain ratios of user-H2O layer-amount to standard value Xh2o = Pwater/Ch2o(Nzc) C ** Calculate ratio of user-defined total ozone to C ** standard-atmosphere column amount C Xo3 = 1.0 IF ( Ozone.GE.0.0 ) Xo3 = Ozone / Co3( Nzc ) C ** Calculate layer amounts of water vapor and ozone DO Lc = 1, Nzc-1 H2oamt( Lc ) = Xh2o * ( Ch2o(Lc+1)-Ch2o(Lc) ) O3amt( Lc ) = Xo3 * ( Co3(Lc+1)-Co3(Lc) ) END DO C ** Re-calculate cumulative water vapor and ozone amounts C Ch2o( 1:Nzc ) = 0.0 Co3( 1:Nzc ) = 0.0 DO Lc = 1, Nzc-1 Ch2o( Lc+1) = Ch2o( Lc ) + H2oamt( Lc ) Co3( Lc+1) = Co3( Lc ) + O3amt( Lc ) END DO C ** Find layers with aerosol and cloud extinction C CALL Mark( Nzc, Zc, Nlaer, Zaert, Zaerb, Laer ) CALL Mark( Nzc, Zc, Nlcld, Zcldt, Zcldb, Lcld ) IF ( Prnt ) THEN PRINT 9000 Lu = 1 DO Lc = 1, Nzc-1 IF ( Lc.EQ.Lprn(Lu) ) THEN PRINT 9020, Lc, Zc(Lc), Pc(Lc), T(Lc), H2oamt(Lc), & Ch2o(Lc), O3amt(Lc), Co3(Lc), Laer(Lc), & Lcld(Lc) Lu = Lu + 1 ELSE IF ( Lc.NE.Lprn(Lu) ) THEN PRINT 9025, Lc, Zc(Lc), Pc(Lc), T(Lc), H2oamt(Lc), & Ch2o(Lc), O3amt(Lc), Co3(Lc), Laer(Lc), & Lcld(Lc) END IF END DO IF ( Nzc.EQ.Lprn(Lu) ) THEN PRINT 9030, Lc, Zc(Nzc), Pc(Nzc), T(Nzc), & Ch2o(Nzc), Co3(Nzc) ELSE IF ( Nzc.NE.Lprn(Lu) ) THEN PRINT 9035, Lc, Zc(Nzc), Pc(Nzc), T(Nzc), & Ch2o(Nzc), Co3(Nzc) END IF END IF 9000 FORMAT( '1', 120('*'), /, 45X, ' ATMOSPHERIC VERTICAL STRUCTURE', & /1X, 120('*'), /, & ' HEIGHT PRESS TEMP LYR H2OAMT', & ' CUM H2OAMT LYR O3AMT CUM O3AMT AEROSOL CLOUD', & /, 8X, & ' (KM) (MB) (K) (G/SQ M) (G/SQ CM)', & ' (ATM-CM) (ATM-CM) M LAYER MODEL' ) 9020 FORMAT ( I3, '**', F9.3, F9.2, F8.2, 1P, 4E11.3, I8, I7 ) 9025 FORMAT ( I3, 2X, F9.3, F9.2, F8.2, 1P, 4E11.3, I8, I7 ) 9030 FORMAT ( I3, '**', F9.3, F9.2, F8.2, 2(11X, 1P, E11.3) ) 9035 FORMAT ( I3, 2X, F9.3, F9.2, F8.2, 2(11X, 1P, E11.3) ) CONTAINS C*********************************************************************** SUBROUTINE Addlev( Nzc, Zc, Nadd, Zadd, Maxlev ) C*********************************************************************** C C Inserts levels -Zadd- into array -Zc- C C INPUT: C C Nadd : number of levels to add to -Z- C Ncz : number of levels C Zadd(L) : L=1 to Nadd; levels (km) to add to -Z- C Zc(L) : L=1 to Nzc; altitude (km) C Maxlev : max number of levels C C OUTPUT : C C Ncz : number of levels C Zc(L) : L=1 to Nzc; altitude (km) including levels -Zadd- C C Inserts levels Zadd into Zc C C .. Implicit None Statement .. IMPLICIT None C .. C .. Scalar Arguments .. INTEGER Maxlev, Nadd, Nzc C .. C .. Array Arguments .. REAL Zadd( * ), Zc( * ) C .. C .. Local Scalars .. LOGICAL Dimerr INTEGER K, L, La, Lc, Nz REAL Thin C .. C .. Local Arrays .. REAL Z( Maxlev ) C .. C .. External Functions .. c LOGICAL Wrtdim c EXTERNAL Wrtdim C .. C .. External Subroutines .. c EXTERNAL Errmsg C .. C .. Data statements .. DATA Thin / 1.E-3 / C .. C .. Executable Statements .. C ** Avoid creation of very thin layers by setting C ** values of Zc to those of Zadd if the difference C ** between them is less than or equal to Thin DO Lc = 1, Nzc DO L = 1, Nadd IF ( ABS(Zadd(L)-Zc(Lc)).LE.Thin ) THEN Zc( Lc ) = Zadd( L ) EXIT END IF END DO END DO C ** Insert new level C Z = 0.0 Nz = Nzc Z( 1:Nz ) = Zc( 1:Nzc ) Lc = 0 DO L = 1, Nz-1 Lc = Lc + 1 Zc( Lc ) = Z( L ) DO La = 1, Nadd IF ( Zadd(La).LT.Z(L) .AND. Zadd(La).GT.Z(L+1) ) THEN Lc = Lc + 1 Zc( Lc ) = Zadd( La ) K = La END IF END DO END DO Nzc = Lc + 1 Dimerr = .FALSE. IF ( Nzc.GT.Maxlev ) Dimerr = Wrtdim( 'Maxlev', Nzc ) IF ( Dimerr ) CALL Errmsg( 'Addlev--Dimension error', .TRUE. ) Zc( Nzc ) = Z( Nz ) END SUBROUTINE Addlev C*********************************************************************** SUBROUTINE Mark( Nz, Z, Nlayer, Ztop, Zbot, Layer ) C*********************************************************************** C C INPUT: C C Nlayer : number of layers C Nz : number of levels (-Z- values) C Z(L ) : L=1 to Nz; altitude (km) C Ztop(La) : La=1 to Nlayer; altitude of top of layer (km) C Zbot(La) : La=1 to Nlayer; altitude of bottom of layer (km) C C OUTPUT: C C Layer(L) : L=1, Nz-1; C .. Implicit None Statement .. IMPLICIT None C .. C .. Scalar Arguments .. INTEGER Nlayer, Nz C .. C .. Array Arguments .. INTEGER Layer( Nz-1 ) REAL Z( Nz ), Zbot( Nlayer ), Ztop( Nlayer ) C .. C .. Local Scalars .. INTEGER L, La REAL Thin C .. C .. Data statements .. DATA Thin / 1.E-3 / C .. C .. Executable Statements .. Layer = 0 DO L = 1, Nz-1 DO La = 1, Nlayer IF ( Z(L).LE.Ztop(La)+Thin .AND. & Z(L+1).GE.Zbot(La)-Thin ) THEN Layer( L ) = La EXIT END IF END DO cc print *, L, Z(L), Z(L+1), Layer( L ) END DO END SUBROUTINE Mark END SUBROUTINE Setatm C*********************************************************************** REAL FUNCTION Averag( F1, F2 ) C*********************************************************************** C C Compute the average value of a function f(x) across a region, C given its values (f1, f2) at the boundaries of the region, C and assuming f(x) varies exponentially unless C (a) f1 is nearly equal to f2, or C (b) f1 or f2 is near zero. C C .. Implicit None Statement .. IMPLICIT None C .. C .. Scalar Arguments .. REAL F1, F2 C .. C .. Local Scalars .. REAL Equal, Ratio, Zero C .. C .. Intrinsic Functions .. INTRINSIC Abs, Log C .. C .. Data statements .. DATA Zero / 1.E-20 / , Equal / 0.01 / C .. C .. Executable Statements .. IF ( F1.GT.Zero .AND. F2.GT.Zero ) THEN Ratio = F2 / F1 IF ( Abs(Ratio-1.0).GE.Equal ) THEN Averag = ( F2-F1 ) / Log( Ratio ) ELSE Averag = 0.5 * ( F1+F2 ) END IF ELSE Averag = 0.5 * ( F1+F2 ) END IF RETURN END C*********************************************************************** REAL FUNCTION Densty( P, T, R ) C*********************************************************************** C C Calculates atmospheric density (kg/cu m) C C P : pressure (mb) C T : temperature (k) C R : mixing ratio (g/g) C RD: gas constant for dry air (mks) C C .. Implicit None Statement .. IMPLICIT None C .. C .. Scalar Arguments .. REAL P, R, T C .. C .. Local Scalars .. REAL Rd C .. C .. External Functions .. c REAL Virtem c EXTERNAL Virtem C .. C .. Data statements .. DATA Rd / 287.05 / C .. C .. Executable Statements .. Densty = 100. * P / ( Rd*Virtem(T,R) ) RETURN END C*********************************************************************** REAL FUNCTION Virtem( T, R ) C*********************************************************************** C C Calculates virtual temperature (K) C C T : temperature (k) C R : mixing ratio (g/g) C .. Implicit None Statement .. IMPLICIT None C .. C .. Scalar Arguments .. REAL R, T C .. C .. Local Scalars .. REAL Eps C .. C .. Data statements .. DATA Eps / 0.622 / C .. C .. Executable Statements .. Virtem = T * ( 1.+R/Eps ) / ( 1.+R ) RETURN END C*********************************************************************** SUBROUTINE Alt( Nlevu, Pu, Z, P, Nz, Maxlev, Zu ) C*********************************************************************** C C Calculates altitude Zu(km) at pressure level -Pu- by C interpolating Z(P) C C INPUT: C C Maxlev : max no. of levels, dimensions of -P- and -Z- C Nlevu : number of levels at which to calculate altitude C Nz : no. of levels in -Z- and -P- C Pu(Lu) : Lu=1 to Nlevu; pressure at which to calculate altitude C Z(L) : L=1 to Nz; altitude C P(L) : L=1 to Nz; pressure C C OUTPUT C C Zu(Lu) : Lu=1 to Nlevu; altitude at pressure -Pu- C C .. Implicit None Statement .. IMPLICIT None C .. C .. Scalar Arguments .. INTEGER Maxlev, Nlevu, Nz C .. C .. Array Arguments .. REAL P( Maxlev ), Pu( Nlevu ), Z( Maxlev ), Zu( Nlevu ) C .. C .. Local Scalars .. INTEGER L, Lu C .. C .. Executable Statements .. DO Lu = 1, Nlevu L = 2 DO WHILE ( P(L).LT.Pu(Lu) .AND. L.LT.Nz ) L = L + 1 END DO Zu( Lu ) = Log10( Pu(Lu)/P(L-1))/Log10(P(L)/P(L-1) ) * & ( Z(L) - Z(L-1) ) + Z( L-1 ) cc print *, Lu, Pu(Lu), Zu(Lu), L-1, L END DO END C*********************************************************************** SUBROUTINE Setcld( Nzc, Zc, Nlcld, Cldtyp,Taucld, Clder, Zcldt, & Zcldb, Lcld, Nwl, Wavlen, Nuvis, Musun, & Maxlev, Maxwvl,Dtauex, Ssalb, Gf ) C*********************************************************************** C C Calculates extinction optical depth, asymmetry factor and C single scattering albedo of cloud in computational layers C C INPUT VARIABLES: C C Lcld(Lc) : Lc=1 to Nzc, model-cloud layer to be used in C computational layer Lc C Maxlev : max no. of vertical levels C Maxwvl : max no. of spectral intervals C Musun : cosine of solar zenith angle C Nlcld : number of cloud layers C Nwl : number of spectral intervals C Nuvis : number of spectral intervals in the uv-visible range C Nzc : number of vertical levels C Cldtyp(L) : L=1 to Nlcld; Cloud type-1:Water, 2:Ice C Taucld(L) : L=1 to Nlcld; optical depth of cloud layer -L- C Clder(L) : L=1 to Nlcld; Droplet effective radius of cloud layer -um- C Wavlen(I) : I=1 to Nwl+1; wavelength (micron); -Wavlen(I)- is the C lower boundary and -Wavlen(I+1)- is the upper boundary C of a spectral interval C Zc(Lc) : Lc=1 to Nzc, height of level -Lc- (km) C Zcldb(L) : L=1 to Nlcld; cloud-base height for cloud layer -L- (km) C Zcldt(L) : L=1 to Nlcld; cloud-top height for cloud layer -L- (km) C C OUTPUT VARIABLES: C C Dtauex(Lc): Lc=1 to Nzc-1, extinction optical depth of computational C layer -Lc- C Gf(Lc) : Lc=1 to Nzc-1, asymmetry factor of computational layer C -Lc- C Ssalb(Lc) : Lc=1 to Nzc-1, single scattering albedo of computational C layer -Lc- C .. Implicit None Statement .. IMPLICIT NONE C .. C .. Scalar Arguments .. INTEGER Maxlev, Maxwvl, Nlcld, Nwl, Nzc,Nuvis REAL Musun C .. C .. Array Arguments .. INTEGER Lcld( * ),Cldtyp( * ) REAL Dtauex( Maxlev-1, Maxwvl ), Gf( Maxlev-1, Maxwvl ), & Ssalb( Maxlev-1, Maxwvl ), Taucld( * ), Wavlen( * ), & Zc( * ), Zcldb( * ), Zcldt( * ), Clder( * ) C .. C .. Local Scalars .. INTEGER I, Lc, N, Nl REAL Tau(Nwl) C .. C .. Local Arrays .. REAL Bncld( Nwl,Nlcld ), Gfcld(Nwl, Nlcld ), Oucld(Nwl, Nlcld ) C .. C .. External Subroutines .. c EXTERNAL Cldwt, Cldic C .. DO N = 1, Nlcld IF ( Cldtyp(N) .EQ. 1) THEN CALL Cldwt( Taucld(N),Clder(N),Oucld(1,N),Gfcld(1,N),Tau ) ELSE CALL Cldic( Taucld(N),Clder(N),Oucld(1,N),Gfcld(1,N),Tau ) END IF DO i=1, Nwl Bncld(i, N ) = Tau( i ) / ( Zcldt(N)-Zcldb(N) ) ENDDO ENDDO Dtauex = 0.0 Ssalb = 0.0 Gf = 0.0 DO 40 Lc = 1, Nzc - 1 IF ( Lcld(Lc).NE.0 ) THEN DO 30 I = 1, Nwl Dtauex( Lc, I ) = Bncld(I, Lcld(Lc) ) * & ( Zc(Lc)-Zc(Lc+1))*1000. !test 120425 c if (dtauex(lc,i)<0) then c print*, '----inside' c print*, lc,i,dtauex(lc,i),bncld(i,lcld(lc)), c & ( Zc(Lc)-Zc(Lc+1))*1000. c endif !test 120425 end Ssalb( Lc, I ) = Oucld(I, Lcld(Lc) ) Gf( Lc, I ) = Gfcld(I, Lcld(Lc) ) 30 CONTINUE END IF 40 CONTINUE !test 120425 c if (any(dtauex<0)) then c print*,'------dtauex' c print*, dtauex c print*, '-----nlcld' c print*, nlcld c print*, '------lcld' c print*, lcld(1:nzc) c print*, '------bncld' c print*, bncld(1:nwl,1:nlcld) c print*, '------zcl-zcl+1' c print*, (zc(1:nzc-1)-zc(2:nzc))*1000. c print*, '------tau' c print*, tau(1:nwl) c print*, '------zcldt' c print*, zcldt(1:nlcld) c print*, '------zcldb' c print*, zcldb(1:nlcld) c print*, '------ssalb' c print*, ssalb(1:Nzc-1,:) c print*, '------dtauex' c print*, dtauex(1:nzc-1,:) c print*, '------gf' c print*, gf(1:nzc-1,:) c endif !test 120425 end END C*********************************************************************** SUBROUTINE Cldwt( Tau, Clder, Oucld, Gfcld, Taucld) C*********************************************************************** C C Calculates Cloud extinction coeffcient, cloud single scattering albedo C and asymmetry factor from the parameterization of Slingo (1989), C and Edwards et al. (1996) for water clouds. C C Input: Tau : visible cloud optical depth C Clder : Cloud droplet effective radius C C Output: Gfcld(i) : i=1,Nwl; asymmetry factor C Oucld(i) : i=1,Nwl; single scattering albedo C Taucld(i) : i=1,Nwl; cloud optical depth C C Internal Variables: C C coef(i,j) : i=1 to 4, j=1 to 6; Coeffcients used to estimate cloud optical properties C .. Implicit None Statement .. IMPLICIT NONE C .. C .. Parameters .. INTEGER Nwl PARAMETER ( Nwl=7 ) C .. C .. Scalar Arguments .. REAL Tau, Clder C .. Array Arguments .. REAL Gfcld(Nwl), Oucld(Nwl), Taucld(Nwl) C .. C .. Local Scalars .. INTEGER I, J, Il C .. C .. Local Arrays .. REAL Coef(4,6) C .. C .. Data statements .. DATA ( Coef(1,I), I=1, 6 ) / -8.737E0, 1.671E-3, 7.465E-8, & 1.114E-1, 8.371E-1, 1.729E3 / DATA ( Coef(2,I), I=1, 6 ) / -1.451E1, 1.772E-3, 6.745E-6, & 9.879E0, 8.144E-1, 2.642E3 / DATA ( Coef(3,I), I=1, 6 ) / -2.576E1, 1.959E-3, 1.278E-3, & 6.149E2, 7.914E-1, 3.701E3 / DATA ( Coef(4,I), I=1, 6 ) / -3.414E1, 2.147E-3, 6.067E-2, & 7.566E3, 8.354E-1, 3.503E3 / C .. C ** Get Cloud optical depth, single scattering albedo and C ** asymmetry factor IF ( clder .LT. 3.0)clder=3.0 IF ( clder .GT. 24.0)clder=24.0 clder=clder*1.0E-6 IF ( Tau.NE.0.0 .AND. Clder .NE.0.0) THEN DO 10 il=1, 4 Taucld( il )=Tau Oucld( il )=1.-(Coef(1,3)+Coef(1,4)*Clder) Gfcld( Il )=Coef(1,5)+Coef(1,6)*Clder 10 CONTINUE DO 20 il=5, Nwl Taucld( il )=Tau*(Coef(il-3,1)+Coef(il-3,2)/Clder)/ & (Coef(1,1)+Coef(1,2)/Clder) Oucld( il )=1.-(Coef(il-3,3)+Coef(il-3,4)*Clder) Gfcld( Il )=Coef(il-3,5)+Coef(il-3,6)*Clder 20 CONTINUE c write(*,*)taucld c write(*,*)oucld c write(*,*)gfcld END IF END C*********************************************************************** SUBROUTINE Cldic( Tau, Clder, Oucld, Gfcld, Taucld) C*********************************************************************** C C Calculates Cloud extinction coeffcient, cloud single scattering albedo C and asymmetry factor from the parameterization of Chou et al.(2002) for ice clouds. C C Input: Tau : visible cloud optical depth C Clder : Cloud droplet effective radius C C Output: Gfcld(i) : i=1,Nwl; asymmetry factor C Oucld(i) : i=1,Nwl; single scattering albedo C Taucld(i) : i=1,Nwl; cloud optical depth C C Internal Variables: C C coef(i,j) : i=1 to 5, j=1 to 6; Coeffcients used to estimate cloud optical properties C .. Implicit None Statement .. IMPLICIT NONE C .. C .. Parameters .. INTEGER Nwl PARAMETER ( Nwl=7 ) C .. C .. Scalar Arguments .. REAL Tau, Clder C .. Array Arguments .. REAL Gfcld(Nwl), Oucld(Nwl), Taucld(Nwl) C .. C .. Local Scalars .. INTEGER I, J, Il C .. C .. Local Arrays .. REAL Coef(5,6) C .. C .. Data statements .. DATA ( Coef(1,I), I=1, 6 ) / 1.37E-7, 7.06E-8, 5.64E-12, & 7.56E-1, 1.08E-3, -4.12E-6 / DATA ( Coef(2,I), I=1, 6 ) / -1.52E-7, 7.38E-8, -3.48E-11, & 7.46E-1, 1.41E-3, -5.74E-6 / DATA ( Coef(3,I), I=1, 6 ) / 1.41E-6, 5.72E-6, -1.22E-9, & 7.25E-1, 1.85E-3, -7.73E-6 / DATA ( Coef(4,I), I=1, 6 ) / 1.12E-3, 5.65E-4, -8.96E-7, & 7.17E-1, 2.28E-3, -8.86E-6 / DATA ( Coef(5,I), I=1, 6 ) / 4.83E-2, 2.74E-3, -9.02E-6, & 7.71E-1, 2.45E-3, -1.00E-5 / C .. C ** Get Cloud optical depth, single scattering albedo and C ** asymmetry factor IF ( Tau.NE.0.0 .AND. Clder .NE.0.0) THEN IF (clder .LT. 10.0) clder=10. IF (clder .GT. 150.0)clder=150. Taucld( 1 )=Tau Oucld( 1 )=1.-(Coef(1,1)+Coef(1,2)*Clder $ +Coef(1,3)*Clder**2) Gfcld( 1 )=Coef(1,4)+Coef(1,5)*Clder+ $ Coef(1,6)*Clder**2 DO 10 il=2, 4 Taucld( il )=Tau Oucld( il )=1.-(Coef(2,1)+Coef(2,2)*Clder+ $ Coef(2,3)*Clder**2) Gfcld( Il )=Coef(2,4)+Coef(2,5)*Clder+ $ Coef(2,6)*Clder**2 10 CONTINUE DO 20 il=5, Nwl Taucld( il )=Tau Oucld( il )=1.-(Coef(il-2,1)+Coef(il-2,2)*Clder+ $ Coef(il-2,3)*Clder**2) Gfcld( il )=Coef(il-2,4)+Coef(il-2,5)*Clder+ $ Coef(il-2,6)*Clder**2 20 CONTINUE END IF c write(*,*)taucld c write(*,*)oucld c write(*,*)gfcld END C*********************************************************************** SUBROUTINE Solar( Iwl, Musun, Rsun, Solcon, Solflx, Flxsun ) C*********************************************************************** C Calculates solar flux for current spectral interval from C MODTRAN-3 C C C INPUT: C Iwl : index of current spectral interval C Musun : cosine of solar zenith angle C Rsun : correction factor for mean solar flux (square of ratio C of mean to actual earth-sun distance) C Solcon : solar constant (watts/m**2) C C OUTPUT: C Flxsun : solar flux normal to the atmosphere (watts/m**2) C Solflx : solar flux normal to beam (watts/m**2) C C LOCAL VARIABLES: C E(I) : I=1 to Nwl, solar flux (watts/m**2) at 1 astronomical C unit (a.u.), integrated between wavelengths C Wvl(I) and Wvl(I+1) C Scltot : total solar constant from integrating solar spectrum C for 0-45 microns C Wvl(I) : I=1 to Nwl+1, tabular wavelength I C .. Implicit None Statement .. IMPLICIT NONE C .. C .. Parameters .. INTEGER Nwl PARAMETER ( Nwl=7 ) C .. C .. Scalar Arguments .. INTEGER Iwl REAL Musun, Rsun, Solcon, Solflx, Flxsun C .. C .. Local Scalars .. REAL Scltot C .. C .. Local Arrays .. REAL E( Nwl ), Wvl( Nwl+1 ) C .. C .. Data statements .. DATA Wvl / 0.2, 0.4, 0.5, 0.6, 0.7, 1.19, 2.38, 4.0 / DATA E / 111.542, 187.369, 186.301, 159.153, 422.938, & 247.112, 40.622 / DATA Scltot / 1366.94 / c DATA Scltot / 1355.04 / C .. C .. Executable Statements .. Solflx = ( Solcon / Scltot ) * Rsun * E( Iwl ) Flxsun = Musun * Solflx END C**********************************************************************C SUBROUTINE Sundat( Iday, Pi, Rsun, Decl ) C**********************************************************************C C C For a given day of the year computes the declination of the C sun, equation of time and the earth-sun distance factor. C C Reference: C Spencer, J. W., 1971: Fourier series representation of the C position of the sun. Search 2 (5), 172 C C I N P U T V A R I A B L E S : C C Iday = day of year (1 on 1 January) C Pi = 3.1415927... C C O U T P U T V A R I A B L E S : C C Decl = declination of sun (in radians) C Rsun = earth-sun distance factor ( square of the ratio of mean C distance to actual one) C Eqtime = equation of time (in hours) C C L O C A L V A R I A B L E: C C Dangle = day angle (in radians) C C C .. Scalar Arguments .. INTEGER Iday REAL Decl, Eqtime, Pi, Rsun C .. C .. Local Scalars .. REAL Dangle C .. C .. Intrinsic Functions .. INTRINSIC Cos, Sin C .. Dangle = 2. * Pi * ( Iday-1 ) / 365. Rsun = 1.00011 + 0.034221 * Cos( Dangle ) + & 0.00128 * Sin( Dangle ) + 0.000719 * Cos( 2. *Dangle ) + & 0.000077 * Sin( 2. *Dangle ) Decl = ( 0.006918-0.399912 *Cos(Dangle)+0.070257 *Sin(Dangle)- & 0.006758 *Cos(2. *Dangle)+0.000907 *Sin(2. *Dangle)- & 0.002697 *Cos(3. *Dangle)+0.00148 *Sin(3. *Dangle) ) Eqtime = ( 0.000075+0.001868 *Cos(Dangle)-0.032077 *Sin(Dangle)- & 0.014615 *Cos(2. *Dangle)-0.04089 *Sin(2. *Dangle) ) * & 12. / Pi RETURN END C**********************************************************************C FUNCTION Sums( Flux, Iw1, Iw2, Nwltot ) C**********************************************************************C C Sums up elements of the NwlTot-element array Flux between C elements Iw1 and Iw2 C INPUT: C Flux(Iv) : Iv=1 to Nwltot; spectral flux C Iw1 : starting index of spectral summation C Iw2 : ending index of spectral summation C Nwltot : total number of spectral intervals in -Flux- C OUTPUT: C Sums : sum of spectral flux between intervals of -Iw1- and -Iw2- C .. Scalar Arguments .. INTEGER Iw1, Iw2, Nwltot, Maxwvl PARAMETER ( Maxwvl=7 ) C .. C .. Array Arguments .. REAL Flux(*) C .. C .. Local Scalars .. INTEGER Iw LOGICAL Inperr C .. C .. Function Return Value .. REAL Sums C .. Inperr = (Iw1 .LT. 1) .OR. (Iw2 .GT. Nwltot) IF (Inperr) CALL Errmsg('Sums-Input error(s)',.TRUE.) Sums = 0.0 DO 10 Iw = Iw1, Iw2 Sums = Sums + Flux(Iw) 10 CONTINUE RETURN END C*********************************************************************** SUBROUTINE Errmsg( Messag, Fatal ) C*********************************************************************** C Print out a warning or error message; abort if error C after making symbolic dump (machine-specific) C .. Scalar Arguments .. CHARACTER Messag * ( * ) LOGICAL Fatal C .. C .. Local Scalars .. LOGICAL Msglim INTEGER Maxmsg, Nummsg C .. C .. Save statement .. SAVE Maxmsg, Nummsg, Msglim C .. C .. Data statements .. DATA Nummsg / 0 / , Maxmsg / 100 / , Msglim / .FALSE. / C .. IF ( Fatal ) THEN WRITE ( *, FMT='(//,2A,//)' ) ' ******* ERROR >>>>>> ', Messag STOP END IF Nummsg = Nummsg + 1 IF ( Msglim ) RETURN IF ( Nummsg.LE.Maxmsg ) THEN WRITE ( *, FMT='(/,2A,/)' ) ' ******* WARNING >>>>>> ', Messag ELSE WRITE ( *, FMT=9000 ) Msglim = .True. END IF RETURN 9000 FORMAT ( /, /, ' >>>>>> TOO MANY WARNING MESSAGES -- ', & 'They will no longer be printed <<<<<<<', /, / ) END C*********************************************************************** LOGICAL FUNCTION Wrtbad( Varnam ) C*********************************************************************** C Write names of erroneous variables and return 'TRUE' C INPUT : VarNam = Name of erroneous variable to be written C ( CHARACTER, any length ) C .. Scalar Arguments .. CHARACTER Varnam * ( * ) C .. C .. Local Scalars .. INTEGER Maxmsg, Nummsg C .. C .. External Subroutines .. c EXTERNAL Errmsg C .. C .. Save statement .. SAVE Nummsg, Maxmsg C .. C .. Data statements .. DATA Nummsg / 0 / , Maxmsg / 50 / C .. Wrtbad = .TRUE. Nummsg = Nummsg + 1 WRITE ( *, FMT='(3A)' ) ' **** Input variable ', Varnam, & ' in error ****' IF ( Nummsg.EQ.Maxmsg ) CALL Errmsg( & 'Too many input errors. Aborting...' & , .TRUE. ) RETURN END C*********************************************************************** LOGICAL FUNCTION Wrtdim( Dimnam, Minval ) C*********************************************************************** C Write name of too-small symbolic dimension and C the value it should be increased to; return 'TRUE' C INPUT : DimNam = Name of symbolic dimension which is too small C ( CHARACTER, any length ) C Minval = Value to which that dimension should be C increased (at least) C .. Scalar Arguments .. CHARACTER Dimnam * ( * ) INTEGER Minval C .. WRITE ( *, FMT='(3A,I7)' ) ' **** Symbolic dimension ', Dimnam, & ' should be increased to at least ', Minval Wrtdim = .TRUE. RETURN END C*********************************************************************** LOGICAL FUNCTION Tstbad( Varnam, Relerr ) C*********************************************************************** C Write name (VarNam) of variable failing self-test and its C percent error from the correct value; return 'FALSE'. C .. Scalar Arguments .. CHARACTER Varnam * ( * ) REAL Relerr C .. Tstbad = .FALSE. WRITE ( *, FMT='(/,3A,1P,E11.2,A)' ) ' Output variable ', Varnam, & ' differed by ', 100. * Relerr, & ' per cent from correct value. Self-test failed.' RETURN END END MODULE ! end of Module_SRB