program extract_arts ! 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 91-level diverse profile dataset from ECMWF ! Owner: ! Frederic Chevallier ! ! History: ! Version Date Comment ! 1 22/11/2006 Frederic Chevallier, LSCE ! Code description: ! Language: Fortran 90. ! Software Standards: "European Standards for Writing and Documenting ! Exchangeable Fortran 90 code". implicit none call do_dataset("t ") call do_dataset("q ") call do_dataset("oz ") call do_dataset("ccol ") call do_dataset("rcol ") end program extract_arts !-------------------------------------------------------------------- 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 DO_DATASET(cvar) ! Do conversion for one of the 5 datasets. implicit none integer, parameter :: nlev = 91 ! 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 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 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 ! write(*,*) "Enter the identification of the sampled variable:" ! write(*,*) "- t (for temperature)" ! write(*,*) "- q (for humidity)" ! write(*,*) "- oz (for ozone)" ! write(*,*) "- ccol (for cloud condensate)" ! write(*,*) "- rcol (for precipitation)" ! read(*,*) cvar ! write(*,*) cvar write(*,*) "Extracting " , cvar, " ..." 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='storage4/home/sbuehler/scalc/2007_10_chevallier_91L/original_data/nwp_saf_'//trim(cvar)//& &'_sampled.atm',form='formatted') open(2,file='storage4/home/sbuehler/scalc/2007_10_chevallier_91L/original_data/nwp_saf_'//trim(cvar)//& &'_sampled.sfc',form='formatted') ! ARTS-1-1 output file: open(3,file='chevallierl91_all_'//trim(cvar)//'.xml',form='formatted') write(3,'(A)') '' write(3,'(A)') ' ' write(3,'(A)') ' ' write(3,'(A)') 'This file contains one matrix for each atmospheric state.' 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)') ' ' i = 0 do read(1,'(833(e13.6,x),i5,3i3,i7,i5)',end=1) & temp(:), &! 1) Temperature [K] (1-91) hum(:), &! 2) Humidity [kg/kg] (92-182) ozo(:), &! 3) Ozone [kg/kg] (183-273) cc(:), &! 4) Cloud Cover [0-1] (274-364) clw(:), &! 5) C Liquid W [kg/kg] (365-455) ciw(:), &! 6) C Ice W [kg/kg] (456-546) rain(:), &! 7) Rain [kg/(m2 *s)] (547-637) snow(:), &! 8) Snow [kg/(m2 *s)] (638-728) w(:), &! 9) Vertical Velocity [Pa/s] (729-819) lnpsurf, &!10) Ln of Surf Pressure in [Pa] (820) z0, &!11) Surface geopotential [m2/s2] (821) tsurf, &!12) Surface Skin Temperature [K] (822) t2m, &!13) 2m Temperature [K] (823) td2m, &!14) 2m Dew point temperature [K] (824) hum2m, &!15) 2m Specific Humidity [kg/kg] (825) u10, &!16) 10m wind speed U component [m/s] (826) v10, &!17) 10m wind speed V component [m/s] (827) stratrsrf, &!18) Stratiform rain at surface [kg/(m2 *s)] (828) convrsrf, &!19) Convective rain at surface [kg/(m2 *s)] (829) snowsurf, &!20) Snow at surface [kg/(m2 *s)] (830) lsm, &!21) Land/sea Mask [0-1] (831) lat, &!22) Latitude [deg] (832) long, &!23) Longitude [deg] (833) year, &!24) Year (834) month, &!25) Month (835) day, &!26) Day (836) step, &!27) Step (837) gpoint, &!28) Grid point [1-843490] (838) ind !29) Index (rank-sorted) (839) read(2,'(26(e13.6,x),i5,3i3,i7,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-843490] (31) ind !%32) Index (rank-sorted) (32) ! compute surface pressure [Pa] psurf = exp(lnpsurf) ! compute surface geometric height [m] elev = z0 / grav ! compute full-level pressure (pap) and half-level pressures (paph) [Pa] ! note that all profile variables above are given at full levels call ec_p91l(psurf,pap,paph) ! Write to ARTS output file: ! 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 specific humidity: call specific_to_vmrh2o(hum2m,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(3,100) vol_clw, vol_ciw, rain(nlev), snow(nlev) !100 format(2ES9.2, 2ES10.2) ! 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) ! Update hydrostatic quantities for next iteration p1 = p2 t1 = t2 h2o1 = h2o2 end do write(3,'(A)') '' i = i + 1 enddo ! Check that there were really 5000 profiles in the dataset: 1 if (i .ne. 5000) then write(*,*) 'Incorrect number of profiles found:',i else write(*,*) 'Correct number of 5000 profiles found.' end if ! Close files: close(1) close(2) write(3,'(A)') '' write(3,'(A)') ' ' close(3) return end SUBROUTINE DO_DATASET ! ---------------------------------------- ! SUBROUTINE EC_P91l(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 2006, EUMETSAT, All Rights Reserved. ! ! Description: ! Computes the 91-level vertical pressure grid ! associated to the input surface pressure ! All pressures are in Pa implicit none integer, parameter :: nlev=91 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.000040, & 3.980832, & 7.387186, & 12.908319, & 21.413612, & 33.952858, & 51.746601, & 76.167656, & 108.715561, & 150.986023, & 204.637451, & 271.356506, & 352.824493, & 450.685791, & 566.519226, & 701.813354, & 857.945801, & 1036.166504, & 1237.585449, & 1463.163940, & 1713.709595, & 1989.874390, & 2292.155518, & 2620.898438, & 2976.302246, & 3358.425781, & 3767.196045, & 4202.416504, & 4663.776367, & 5150.859863, & 5663.156250, & 6199.839355, & 6759.727051, & 7341.469727, & 7942.926270, & 8564.624023, & 9208.305664, & 9873.560547, & 10558.881836, & 11262.484375, & 11982.662109, & 12713.897461, & 13453.225586, & 14192.009766, & 14922.685547, & 15638.053711, & 16329.560547, & 16990.623047, & 17613.281250, & 18191.029297, & 18716.968750, & 19184.544922, & 19587.513672, & 19919.796875, & 20175.394531, & 20348.916016, & 20434.158203, & 20426.218750, & 20319.011719, & 20107.031250, & 19785.357422, & 19348.775391, & 18798.822266, & 18141.296875, & 17385.595703, & 16544.585938, & 15633.566406, & 14665.645508, & 13653.219727, & 12608.383789, & 11543.166992, & 10471.310547, & 9405.222656, & 8356.252930, & 7335.164551, & 6353.920898, & 5422.802734, & 4550.215820, & 3743.464355, & 3010.146973, & 2356.202637, & 1784.854614, & 1297.656128, & 895.193542, & 576.314148, & 336.772369, & 162.043427, & 54.208336, & 6.575628, & 0.003160, & 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.000014, & 0.000055, & 0.000131, & 0.000279, & 0.000548, & 0.001000, & 0.001701, & 0.002765, & 0.004267, & 0.006322, & 0.009035, & 0.012508, & 0.016860, & 0.022189, & 0.028610, & 0.036227, & 0.045146, & 0.055474, & 0.067316, & 0.080777, & 0.095964, & 0.112979, & 0.131935, & 0.152934, & 0.176091, & 0.201520, & 0.229315, & 0.259554, & 0.291993, & 0.326329, & 0.362203, & 0.399205, & 0.436906, & 0.475016, & 0.513280, & 0.551458, & 0.589317, & 0.626559, & 0.662934, & 0.698224, & 0.732224, & 0.764679, & 0.795385, & 0.824185, & 0.850950, & 0.875518, & 0.897767, & 0.917651, & 0.935157, & 0.950274, & 0.963007, & 0.973466, & 0.982238, & 0.989153, & 0.994204, & 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_p91l