program extract_arts_l137 ! This program dumps the Eresmaa/McNally 137-level profile dataset to ARTS-XML format. ! The source code is combined by using extract_arts_all.f90 (by SAB and DK) from the provided Chevallier datasets and readsaf137.f90 ! ! John Vinzenz Mrziglod 2015-09-01 ! Modified by Daniel Kreyling 12/12/2010 ! ! All parameters will be written to one single file per data set. ! The files are denoted by "_all_". E.g.: chevallierl91_all_t.xml ! ! Simple program to dump the variables we want from the Chevallier ! data set to an ARTS data file. Create a modified copy of this if ! you want additional or different variables. ! ! This is the version that extracts also cloud parameters. ! ! Original header from Frederic below. I had to change real to real*8 ! to make the reading routine work. ! ! SAB 2007-07-19 ! ! Description: ! ! Read 137-level diverse profile dataset from ECMWF ! Owner: ! Frederic Chevallier ! ! History: ! Version Date Comment ! 1 22/11/2006 Frederic Chevallier, LSCE ! 1.1 02/09/2014 Update to 137 levels, Reima Eresmaa, ECMWF ! ! Code description: ! Language: Fortran 90. ! Software Standards: "European Standards for Writing and Documenting ! Exchangeable Fortran 90 code". implicit none call readsaf137("t ", 0) call readsaf137("q ", 1) call readsaf137("oz ", 2) call readsaf137("ccol ", 3) call readsaf137("rcol ", 4) end program extract_arts_l137 !-------------------------------------------------------------------- SUBROUTINE dewpoint_temp_to_vmrh2o(Td, p, x) ! Convert the dew point temperature and air pressure to volume mixing ratio x. ! FIXME: I am not sure about this formula and its result (JVM) implicit none real*8 :: Td ! input real*8 :: p ! input real*8 :: x ! output ! Coefficients for Ew: real*8 :: a = -6096.9385 real*8 :: b = 21.2409642 real*8 :: c = -2.711193e-2 real*8 :: d = 1.673952e-5 real*8 :: e = 2.433502 real*8 :: ew ew = exp( a/Td + b + c*Td + d*Td**2 + e*log(Td) ) ! I worked out this formula after 3.65 of W&H, second edition: x = ew / p return end SUBROUTINE dewpoint_temp_to_vmrh2o SUBROUTINE specific_to_vmrh2o(q,x) ! Convert specific humidity q [kg/kg] to volume mixing ratio x. implicit none ! molecular mass of water vapor relative to dry air, from Wallace&Hobbs: real*8, parameter :: eps=0.622 real*8 :: q ! input real*8 :: x ! output ! I worked out this formula after 3.59 of W&H, second edition: x = q / ( q + eps*(1-q) ) ! It is tempting to divide out q, but then the formula is ! numerically unstable for small q (division by zero). ! ! There is also a simplified formula which could be used ! approximately for small q values: ! x = q / ( eps*(1-q) ) ! The difference is that in the correct formula the q ! shows up in the denominator again. return end SUBROUTINE specific_to_vmrh2o !-------------------------------------------------------------------- SUBROUTINE specific_to_vmro3(q,x) ! Convert specific ozone q [kg/kg] to volume mixing ratio x. ! This is a clone from the similar subroutine for h2o. implicit none ! atomic mass of oxygen: mo3 = 15.9994 ! (Ulshöfer und Hornschuh Formelsammlung, 3. Auflage) ! effectiv molecular weight of air: mair = 28.97 (Wallace and Hobbs, 2nd edition) ! ! Molecular weight of O3 relative to (dry) air: ! eps = 3*mo3 / mair real*8, parameter :: eps=1.657 real*8 :: q ! input real*8 :: x ! output ! I worked out this formula after 3.59 of W&H, second edition: x = q / ( q + eps*(1-q) ) ! It is tempting to divide out q, but then the formula is ! numerically unstable for small q (division by zero). ! ! There is also a simplified formula which could be used ! approximately for small q values: ! x = q / ( eps*(1-q) ) ! The difference is that in the correct formula the q ! shows up in the denominator again. return end SUBROUTINE specific_to_vmro3 !-------------------------------------------------------------------- SUBROUTINE condensate_kg_per_kg_to_kg_per_meter3(in, p, T, out) ! Convert condensate in kg/kg to kg/m^3 ! in : condendaste in kg/kg ! p : pressure in Pa ! T : temperature in K ! out : condensate in kg/m^3 ! ! I calculate the volume of 1 kg of air to do the conversion. For ! this, we assume dry air, i.e., I neglect the different mass of ! the H2O that will be also present. I make this simple choice, ! because I am not sure which definition is used inside the ECMWF ! model. implicit none real*8 :: in ! input real*8 :: p ! input real*8 :: T ! input real*8 :: out ! output ! Gas constant of dry air, according to W&H, second edition: real*8 :: Rd = 287 ! Calculate volume of 1kg of dry air, according to ideal gas law: real*8 :: V V = Rd*T/p ! We have to divide by that to go from kg/kg to kg/m^3: out = in / V return end SUBROUTINE condensate_kg_per_kg_to_kg_per_meter3 !-------------------------------------------------------------------- SUBROUTINE delta_z(p1,t1,h2o1,p2,t2,h2o2,dz) ! dz : Vertical thickness (in m) between two (p,t,h2o) gridpoints ! p1,p2 : Pressure [Pa] ! t1,t2 : Temperature [K] ! h2o1,h2o2 : Humidity [VMR] real*8 :: p1,t1,h2o1,p2,t2,h2o2 ! input real*8 :: dz ! output ! molecular mass of water vapor relative to dry air, from Wallace&Hobbs: real*8, parameter :: eps=0.622 ! Gravitational constant (in m/s2) taken from Chevalliers f90 routine real*8, parameter :: grav=9.80665 ! Gas constant of dry air in J/(K kg) from Wallace&Hobbs: real*8, parameter :: Rd=287 ! Virtual temperatures [K]: real*8 :: Tv1, Tv2, Tv ! Scale height [m]: real*8 :: H1,H2,H ! Calculate virtual mean temperature of the layer. ! Virtual temperature from W&H 2nd edition: ! Tv = T / ( 1 - e/p * (1-eps)) ! Where e and p are humidity partial pressure and total pressure, ! respectively. (Thus e/p is the H2O VMR.) Tv1 = t1 / (1 - h2o1*(1-eps)) Tv2 = t2 / (1 - h2o2*(1-eps)) Tv = 0.5 * (Tv1+Tv2) ! Calculate scale height. ! W&H define it as H = Rd/grav * T H1 = Rd/grav * t1 H2 = Rd/grav * t2 H = 0.5 * (H1+H2) dz = H * log(p1/p2) return end SUBROUTINE delta_z !-------------------------------------------------------------------- SUBROUTINE wind_speed_and_direction(vx,vy,wspeed,wdir) ! wspeed : wind speed [m/s] ! wdir : wind direction [-180 to 180 degrees] ! Zero degree is north. Physical flow direction. ! This means, the direction is to where the wind blows. ! ! vx, vy : 2d vector components of the wind implicit none real*8 :: vx, vy ! input real*8 :: wspeed, wdir ! output wspeed = sqrt(vx**2 + vy**2) ! calculate the wind direction wdir = -atan2(vx,vy)*180/3.141592653589793; return end SUBROUTINE wind_speed_and_direction !-------------------------------------------------------------------- SUBROUTINE readsaf137(cvar, findex) implicit none integer, parameter :: nlev = 137 ! number of vertical levels ! Gravitational constant (in m/s2) taken from Chevalliers f90 routine real*8, parameter :: grav=9.80665 real*8, dimension(nlev) :: temp, hum, ozo, cc, clw, ciw real*8, dimension(nlev) :: rain, snow, w, pap real*8, dimension(nlev+1) :: paph real*8 :: lnpsurf, psurf, z0, tsurf, t2m, td2m, hum2m, u10, v10, wspeed, wdirection real*8 :: stratrsrf, convrsrf, snowsurf, elev real*8 :: lsm, lat, long real*8 :: alb, sr, tvl, tvh, cvl, cvh, seaice real*8 :: asn, rsn, tsn, sd, stl1, stl2, stl3, stl4 real*8 :: swvl1, swvl2, swvl3, swvl4, istl1, istl2, istl3, istl4 integer :: year, month, day, step, gpoint, ind, i, j, findex character(len=10) :: cvar ! VMR parameters of H2O and O3: real*8 :: vmrh2o ! H2O Volume mixing ratio real*8 :: vmro3 ! O3 Volume mixing ratio ! Condensate parameters per volume: real*8 :: vol_clw real*8 :: vol_ciw ! For vertical thickness calculation: real*8 :: p1,t1,h2o1,p2,t2,h2o2,dz,zprof if (trim(cvar) .ne. 't' .and. trim(cvar) .ne. 'q' .and. trim(cvar) .ne. 'oz' .and. & trim(cvar) .ne. 'ccol' .and. trim(cvar) .ne. 'rcol') then write(*,*) 'Wrong variable, start again' stop endif open(1,file='profiles137/nwp_saf_'//trim(cvar)//'_sampled.atm',form='formatted') open(2,file='profiles137/nwp_saf_'//trim(cvar)//'_sampled.sfc',form='formatted') ! ARTS-1-1 output file atmosphere data open(3,file='eresmaal137_all_'//trim(cvar)//'.xml.temp',form='formatted') write(3,'(A)') '' write(3,'(A)') ' ' write(3,'(A)') ' ' write(3,'(A)') 'This file contains one matrix for each atmospheric profile.' write(3,'(A)') 'Each matrix row corresponds to one pressure level. The ' write(3,'(A)') 'meaning of the columns is: ' ! write(3,'(A)') 'p[Pa] T[K] z[m] H2O_VMR[1] O3_VMR[1] CLW[kg/m^3] CIW[kg/m^3] Rain[kg/(m2*s)] Snow[kg/(m2*s)]' write(3,'(A)') 'p[Pa] T[K] z[m] CLW[kg/m^3] CIW[kg/m^3] Rain[kg/(m2*s)] Snow[kg/(m2*s)] H2O_VMR[1] O3_VMR[1]' ! write(3,'(A)') 'CLW[kg/m^3] CIW[kg/m^3] Rain[kg/(m2*s)] Snow[kg/(m2*s)]' write(3,'(A)') ' ' write(3,'(A)') ' ' ! ARTS-1-1 output file surface variables and meta data open(4,file='eresmaal137_all_'//trim(cvar)//'_surface.xml',form='formatted') write(4,'(A)') '' write(4,'(A)') ' ' write(4,'(A)') ' ' write(4,'(A)') 'This file contains one row for each atmospheric profile. ' write(4,'(A)') 'The meaning of the columns is:' write(4,'(A)') '1) profile index [1-25000]'//NEW_LINE('A')//'& &2) Ln of Surf Pressure in [Pa]'//NEW_LINE('A')//'& &3) Surface geopotential [m2/s2] '//NEW_LINE('A')//'& &4) Surface Skin Temperature [K]'//NEW_LINE('A')//'& &5) 2m Temperature [K]'//NEW_LINE('A')//'& &6) 2m Dew point temperature [K]'//NEW_LINE('A')//'& &7) 10m wind speed [m/s]'//NEW_LINE('A')//'& &8) 10m wind direction [0-360 DEG]'//NEW_LINE('A')//'& &9) Stratiform precipitation at surface [m]'//NEW_LINE('A')//'& &10) Convective precipitation at surface [m]'//NEW_LINE('A')//'& &11) Snowfall at surface [m] (water equival.)'//NEW_LINE('A')//'& &12) Land/sea Mask [0-1]'//NEW_LINE('A')//'& &13) Latitude [deg]'//NEW_LINE('A')//'& &14) Longitude [deg]'//NEW_LINE('A')//'& &15) Year'//NEW_LINE('A')//'& &16) Month '//NEW_LINE('A')//'& &17) Day '//NEW_LINE('A')//'& &18) Albedo [0-1] '//NEW_LINE('A')//'& &19) Roughness [m] '//NEW_LINE('A')//'& &20) Low vegetation type [0-20] '//NEW_LINE('A')//'& &21) High vegetation type [0-19] '//NEW_LINE('A')//'& &22) Low vegetation fractional cover [0-1]'//NEW_LINE('A')//'& &23) High vegetation fractional cover [0-1]'//NEW_LINE('A')//'& &24) SeaIce cover [0-1] '//NEW_LINE('A')//'& &25) Snow albedo [0-1] '//NEW_LINE('A')//'& &26) Snow density [kg/m3] '//NEW_LINE('A')//'& &27) Snow temperature [K] '//NEW_LINE('A')//'& &28) Snow depth [m] '//NEW_LINE('A')//'& &29) Soil Layer 1 temperature [K] '//NEW_LINE('A')//'& &30) Soil Layer 2 temperature [K] '//NEW_LINE('A')//'& &31) Soil Layer 3 temperature [K] '//NEW_LINE('A')//'& &32) Soil Layer 4 temperature [K] '//NEW_LINE('A')//'& &33) Volumetric Soil Water Layer 1 [m3/m3] '//NEW_LINE('A')//'& &34) Volumetric Soil Water Layer 2 [m3/m3] '//NEW_LINE('A')//'& &35) Volumetric Soil Water Layer 3 [m3/m3] '//NEW_LINE('A')//'& &36) Volumetric Soil Water Layer 4 [m3/m3] '//NEW_LINE('A')//'& &37) Ice temperature Layer 1 [K] '//NEW_LINE('A')//'& &38) Ice temperature Layer 2 [K] '//NEW_LINE('A')//'& &39) Ice temperature Layer 3 [K] '//NEW_LINE('A')//'& &40) Ice temperature Layer 4 [K]' write(4,'(A)') ' ' !write(4,'(A)') '' !write(4,'(A)') 'pindex[1-25000] psurf[Pa] Tsurf[K] wspeed[m/s] wdir[0-360°] land-sea[0-1] sea-ice[0-1] snow-depth[m]' !write(4,'(A)') '' write(4,'(A)') '' ! ARTS-1-1 output file fractional cloud cover open(5,file='eresmaal137_all_'//trim(cvar)//'_extra.xml.temp',form='formatted') write(5,'(A)') '' write(5,'(A)') ' ' write(5,'(A)') ' ' write(5,'(A)') 'This file contains one matrix for each atmospheric profile. ' write(5,'(A)') 'The meaning of column is:' write(5,'(A)') 'p[Pa] cc[0-1] w[m/s]' write(5,'(A)') ' ' write(5,'(A)') ' ' i = 0 do ! read the atmospheric profile read(1,'(1246(e13.6,x),i5,3i3,i7,i5)',end=1) & temp(:), &! 1) Temperature [K] (1-137) hum(:), &! 2) Humidity [kg/kg] (138-274) ozo(:), &! 3) Ozone [kg/kg] (275-411) cc(:), &! 4) Cloud Cover [0-1] (412-548) clw(:), &! 5) C Liquid W [kg/kg] (549-685) ciw(:), &! 6) C Ice W [kg/kg] (686-822) rain(:), &! 7) Rain [kg/(m2 *s)] (823-959) snow(:), &! 8) Snow [kg/(m2 *s)] (960-1096) w(:), &! 9) Vertical Velocity [Pa/s] (1097-1233) lnpsurf, &!10) Ln of Surf Pressure in [Pa] (1234) z0, &!11) Surface geopotential [m2/s2] (1235) tsurf, &!12) Surface Skin Temperature [K] (1236) t2m, &!13) 2m Temperature [K] (1237) td2m, &!14) 2m Dew point temperature [K] (1238) u10, &!15) 10m wind speed U component [m/s] (1239) v10, &!16) 10m wind speed V component [m/s] (1240) stratrsrf, &!17) Stratiform precipitation at surface [m] (1241) convrsrf, &!18) Convective precipitation at surface [m] (1242) snowsurf, &!19) Snowfall at surface [m] (water equival.) (1243) lsm, &!20) Land/sea Mask [0-1] (1244) lat, &!21) Latitude [deg] (1245) long, &!22) Longitude [deg] (1246) year, &!23) Year (1247) month, &!24) Month (1248) day, &!25) Day (1249) step, &!26) Step (1250) gpoint, &!27) Grid point [1-2140702] (1251) ind !28) Index (rank-sorted) (1252) read(2,'(26(e13.6,x),i5,3i3,i8,i5)',end=1) & alb, &! 1) Albedo [0-1] (1) sr, &! 2) Roughness [m] (2) tvl, &! 3) Low vegetation type [0-20] (3) tvh, &! 4) High vegetation type [0-19] (4) cvl, &! 5) Low vegetation fractional cover [0-1] (5) cvh, &! 6) High vegetation fractional cover [0-1] (6) seaice, &! 7) SeaIce cover [0-1] (7) asn, &! 8) Snow albedo [0-1] (8) rsn, &! 9) Snow density [kg/m3] (9) tsn, &!10) Snow temperature [K] (10) sd, &!11) Snow depth [m] (11) stl1, &!12) Soil Layer 1 temperature [K] (12) stl2, &!13) Soil Layer 2 temperature [K] (13) stl3, &!14) Soil Layer 3 temperature [K] (14) stl4, &!15) Soil Layer 4 temperature [K] (15) swvl1, &!16) Volumetric Soil Water Layer 1 [m3/m3] (16) swvl2, &!17) Volumetric Soil Water Layer 2 [m3/m3] (17) swvl3, &!18) Volumetric Soil Water Layer 3 [m3/m3] (18) swvl4, &!19) Volumetric Soil Water Layer 4 [m3/m3] (19) istl1, &!20) Ice temperature Layer 1 [K] (20) istl2, &!21) Ice temperature Layer 2 [K] (21) istl3, &!22) Ice temperature Layer 3 [K] (22) istl4, &!23) Ice temperature Layer 4 [K] (23) ! The following variables also appear in .atm file lsm, &!%24) Land/sea Mask [0-1] (24) lat, &!%25) Latitude [deg] (25) long, &!%26) Longitude[deg] (26) year, &!%27) Year of start of forecast (27) month, &!%28) Month of start of forecast (28) day, &!%29) Day of start of forecast (29) step, &!%30) Forecast step [hours] (30) gpoint, &!%31) Grid point [1-2140702] (31) ind !%32) Index (rank-sorted) (32) ! compute surface pressure [Pa] psurf = exp(lnpsurf) ! compute wind speed and direction call wind_speed_and_direction(u10, v10, wspeed, wdirection) ! Write to ARTS output file (surface variables and meta data): write(4,'(I5,13(e13.6,x),i5,2i3,23(e13.6,x))') & findex*5000+i+1, &!1) file index lnpsurf, &!2) Ln of Surf Pressure in [Pa] z0, &!3) Surface geopotential [m2/s2] tsurf, &!4) Surface Skin Temperature [K] t2m, &!5) 2m Temperature [K] td2m, &!6) 2m Dew point temperature [K] wspeed, &!7) 10m wind speed [m/s] wdirection, &!8) 10m wind direction [0-360 DEG] stratrsrf, &!9) Stratiform precipitation at surface [m] convrsrf, &!10) Convective precipitation at surface [m] snowsurf, &!11) Snowfall at surface [m] (water equival.) lsm, &!12) Land/sea Mask [0-1] lat, &!13) Latitude [deg] long, &!14) Longitude [deg] year, &!15) Year month, &!16) Month day, &!17) Day alb, &!18) Albedo [0-1] sr, &!19) Roughness [m] tvl, &!20) Low vegetation type [0-20] tvh, &!21) High vegetation type [0-19] cvl, &!22) Low vegetation fractional cover [0-1] cvh, &!23) High vegetation fractional cover [0-1] seaice, &!24) SeaIce cover [0-1] asn, &!25) Snow albedo [0-1] rsn, &!26) Snow density [kg/m3] tsn, &!27) Snow temperature [K] sd, &!28) Snow depth [m] stl1, &!29) Soil Layer 1 temperature [K] stl2, &!30) Soil Layer 2 temperature [K] stl3, &!31) Soil Layer 3 temperature [K] stl4, &!32) Soil Layer 4 temperature [K] swvl1, &!33) Volumetric Soil Water Layer 1 [m3/m3] swvl2, &!34) Volumetric Soil Water Layer 2 [m3/m3] swvl3, &!35) Volumetric Soil Water Layer 3 [m3/m3] swvl4, &!36) Volumetric Soil Water Layer 4 [m3/m3] istl1, &!37) Ice temperature Layer 1 [K] istl2, &!38) Ice temperature Layer 2 [K] istl3, &!39) Ice temperature Layer 3 [K] istl4 !40) Ice temperature Layer 4 [K] ! write(4,101) findex*5000+i+1, lnpsurf, tsurf, wspeed, wdirection, lsm, seaice, sd, lat, long !101 format(I5, ES11.4, F7.2, F7.2, F7.2, 2F10.6, F7.2, 2F10.3) ! compute surface geometric height [m] elev = z0 / 9.80665 ! compute full-level pressure (pap) and half-level pressures (paph) [Pa] ! note that all profile variables above are given at full levels call ec_p137l(psurf,pap,paph) ! Write to ARTS output file (atmospheric profile): ! We first output the quantities for the surface. Profile ! quantities given by Chevallier are layer means, not point ! values! ! So, we add another point at the bottom, then simply ! interpret the profile values as if they were point values. ! For the surface, we take the 2m temperature and humidity. We ! do not have a 2m O3 value, so we simply copy the value from ! the lowest layer. ! Calculate vmrh2o from the dew point temperature and air pressure call dewpoint_temp_to_vmrh2o(td2m,psurf,vmrh2o) ! Calculate vmro3 from specific ozone: call specific_to_vmro3(ozo(nlev),vmro3) ! The 4 condensate parameters are also copied from the lowest ! layer. However, we use psurf and t2m for the conversion from ! kg/kg to kg/m^3. ! Calculate cloud liquid water content: call condensate_kg_per_kg_to_kg_per_meter3(clw(nlev), psurf, t2m, vol_clw) ! Calculate cloud ice water content: call condensate_kg_per_kg_to_kg_per_meter3(ciw(nlev), psurf, t2m, vol_ciw) write(3,'(A)') '' ! write(3,100) psurf, t2m, elev, vmrh2o, vmro3, & ! & vol_clw, vol_ciw, rain(nlev), snow(nlev) !100 format(ES11.4, F7.2, F9.2, 2ES9.2, 2ES9.2, 2ES10.2) write(3,100) psurf, t2m, elev, vol_clw, vol_ciw, rain(nlev), snow(nlev), vmrh2o, vmro3 100 format(ES11.4, F7.2, F9.2, 2ES9.2, 2ES10.2, 2ES9.2) ! Write to ARTS output file (fractional cloud cover): write(5,'(A)') '' ! Initialize quantities for hydrostatic z profile: zprof = elev p1 = psurf t1 = t2m h2o1 = vmrh2o do j=nlev,1,-1 ! Calculate vmrh2o from specific humidity: call specific_to_vmrh2o(hum(j),vmrh2o) ! Calculate vmro3 from specific ozone: call specific_to_vmro3(ozo(j),vmro3) ! Quantities for hydrostatic z profile: p2 = pap(j) t2 = temp(j) h2o2 = vmrh2o ! Calculate thickness between this point and the last: call delta_z(p1,t1,h2o1,p2,t2,h2o2,dz) zprof = zprof + dz ! Calculate cloud liquid water content: call condensate_kg_per_kg_to_kg_per_meter3(clw(j), pap(j), temp(j), vol_clw) ! Calculate cloud ice water content: call condensate_kg_per_kg_to_kg_per_meter3(ciw(j), pap(j), temp(j), vol_ciw) ! write(3,100) pap(j), temp(j), zprof, vmrh2o, vmro3, & ! & vol_clw, vol_ciw, rain(j), snow(j) write(3,100) pap(j), temp(j), zprof, vol_clw, vol_ciw, rain(j), snow(j), vmrh2o, vmro3 ! write(3,100) vol_clw, vol_ciw, rain(j), snow(j) ! output fractional cloud cover write(5,'(ES11.4, F10.6, F10.6)') pap(j), cc(j), w(j) ! Update hydrostatic quantities for next iteration p1 = p2 t1 = t2 h2o1 = h2o2 end do write(3,'(A)') '' write(5,'(A)') '' i = i + 1 enddo 1 write(*,*) 'Number of profiles found in the files:',i ! Close files: close(1) close(2) write(3,'(A)') '' write(3,'(A)') ' ' close(3) write(4,'(A)') '' write(4,'(A)') ' ' close(4) write(5,'(A)') '' write(5,'(A)') ' ' close(5) end SUBROUTINE readsaf137 ! ---------------------------------------- ! SUBROUTINE EC_P137l(spres,pap,paph) ! ! This software was developed within the context of ! the EUMETSAT Satellite Application Facility on ! Numerical Weather Prediction (NWP SAF), under the ! Cooperation Agreement dated 25 November 1998, between ! EUMETSAT and the Met Office, UK, by one or more partners ! within the NWP SAF. The partners in the NWP SAF are ! the Met Office, ECMWF, KNMI and MeteoFrance. ! ! Copyright 2014, EUMETSAT, All Rights Reserved. ! ! Description: ! Computes the 137-level vertical pressure grid ! associated to the input surface pressure ! All pressures are in Pa implicit none integer, parameter :: nlev=137 integer :: jk real*8 :: spres real*8 :: aam(nlev+1), bbm(nlev+1) real*8 :: pap(nlev), paph(nlev+1) data aam / & 0.000000, & 2.000365, & 3.102241, & 4.666084, & 6.827977, & 9.746966, & 13.605424, & 18.608931, & 24.985718, & 32.985710, & 42.879242, & 54.955463, & 69.520576, & 86.895882, & 107.415741, & 131.425507, & 159.279404, & 191.338562, & 227.968948, & 269.539581, & 316.420746, & 368.982361, & 427.592499, & 492.616028, & 564.413452, & 643.339905, & 729.744141, & 823.967834, & 926.344910, & 1037.201172, & 1156.853638, & 1285.610352, & 1423.770142, & 1571.622925, & 1729.448975, & 1897.519287, & 2076.095947, & 2265.431641, & 2465.770508, & 2677.348145, & 2900.391357, & 3135.119385, & 3381.743652, & 3640.468262, & 3911.490479, & 4194.930664, & 4490.817383, & 4799.149414, & 5119.895020, & 5452.990723, & 5798.344727, & 6156.074219, & 6526.946777, & 6911.870605, & 7311.869141, & 7727.412109, & 8159.354004, & 8608.525391, & 9076.400391, & 9562.682617, & 10065.978516, & 10584.631836, & 11116.662109, & 11660.067383, & 12211.547852, & 12766.873047, & 13324.668945, & 13881.331055, & 14432.139648, & 14975.615234, & 15508.256836, & 16026.115234, & 16527.322266, & 17008.789062, & 17467.613281, & 17901.621094, & 18308.433594, & 18685.718750, & 19031.289062, & 19343.511719, & 19620.042969, & 19859.390625, & 20059.931641, & 20219.664062, & 20337.863281, & 20412.308594, & 20442.078125, & 20425.718750, & 20361.816406, & 20249.511719, & 20087.085938, & 19874.025391, & 19608.572266, & 19290.226562, & 18917.460938, & 18489.707031, & 18006.925781, & 17471.839844, & 16888.687500, & 16262.046875, & 15596.695312, & 14898.453125, & 14173.324219, & 13427.769531, & 12668.257812, & 11901.339844, & 11133.304688, & 10370.175781, & 9617.515625, & 8880.453125, & 8163.375000, & 7470.343750, & 6804.421875, & 6168.531250, & 5564.382812, & 4993.796875, & 4457.375000, & 3955.960938, & 3489.234375, & 3057.265625, & 2659.140625, & 2294.242188, & 1961.500000, & 1659.476562, & 1387.546875, & 1143.250000, & 926.507812, & 734.992188, & 568.062500, & 424.414062, & 302.476562, & 202.484375, & 122.101562, & 62.781250, & 22.835938, & 3.757813, & 0.000000, & 0.000000 / data bbm / & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000000, & 0.000007, & 0.000024, & 0.000059, & 0.000112, & 0.000199, & 0.000340, & 0.000562, & 0.000890, & 0.001353, & 0.001992, & 0.002857, & 0.003971, & 0.005378, & 0.007133, & 0.009261, & 0.011806, & 0.014816, & 0.018318, & 0.022355, & 0.026964, & 0.032176, & 0.038026, & 0.044548, & 0.051773, & 0.059728, & 0.068448, & 0.077958, & 0.088286, & 0.099462, & 0.111505, & 0.124448, & 0.138313, & 0.153125, & 0.168910, & 0.185689, & 0.203491, & 0.222333, & 0.242244, & 0.263242, & 0.285354, & 0.308598, & 0.332939, & 0.358254, & 0.384363, & 0.411125, & 0.438391, & 0.466003, & 0.493800, & 0.521619, & 0.549301, & 0.576692, & 0.603648, & 0.630036, & 0.655736, & 0.680643, & 0.704669, & 0.727739, & 0.749797, & 0.770798, & 0.790717, & 0.809536, & 0.827256, & 0.843881, & 0.859432, & 0.873929, & 0.887408, & 0.899900, & 0.911448, & 0.922096, & 0.931881, & 0.940860, & 0.949064, & 0.956550, & 0.963352, & 0.969513, & 0.975078, & 0.980072, & 0.984542, & 0.988500, & 0.991984, & 0.995003, & 0.997630, & 1.000000 / do jk=1,nlev+1 paph(jk)=aam(jk)+bbm(jk)*spres end do do jk=1,nlev pap(jk)=0.5*(paph(jk)+paph(jk+1)) end do return end subroutine ec_p137l