Home

The TOMCAT/SLIMCAT Off-Line 3D CTM User's Manual Martyn

image

Contents

1. Ps _ f V v 2 Integrating the continuity equation 1 from the top of the atmosphere to level gives the equation for the vertical velocity In pressure coordinates this is f V v PP Vh Vp 3 In hybrid vertical coordinates the vertical velocity is Op aay OP op 4 AE 4 The divergence is defined as 1 1 U OV DT 5 The relative vorticity is defined as 1 1 OV OU EE 6 atar 6 The streamfunction 4 and the velocity potential x are given by Sk OW dx vat ayes Ah 7 _ 1f 2 OX C Vx 9 D V x 10 4 Model Grid 4 1 Horizontal grid The model grid is variable and determined usually by computational considerations The longitudinal spacing is regular although the latitudinal spacing can be irregular 4 2 Vertical Grid Many GCMs are formulated using a hybrid o p vertical coordinate In this scheme the model levels vary from pure terrain following levels near the surface to pure pressure levels at higher altitudes typically above 100 hPa The pressure p of a model half level interface is then given by Pret Aps t Bpo 11 Consequently many datasets of meteorological analyses are also produced on these hybrid o p levels Appendix 1 defines the symbols used in this report In TOMCAT SLIMCAT the vertical coordinate of the model levels can be defined in two ways i hybrid o p levels or ii hybrid o levels For hybrid p levels the definit
2. cloud base cloud base v V pw dz J Pw V V V Pu v dz surface surface Tf this integral is positive the grid column is convectively unstable and deep convection occurs However if there is a divergence of moisture a further check is made to determine whether the column is unstable with respect to shallow convection Shallow convection depends less on large scale moisture convergence than on evaporation of moisture from the surface If addition of this local surface evaporative flux to the integral above results in a positive water vapour flux at cloud base then shallow convection occurs within the column If there is still a divergence of moisture below cloud base the model column is set as cloudless To maintain moisture content during convection a moisture balance is imposed below cloud base The updraft mass flux of air through cloud base is therefore calculated for shallow and deep convection as cloud base S v Vpydz surface evaporation M s surface dp de at cloud base and the updraft cloud base mass flux of tracer Trbese as Peer Mice con The updraft mass flux at cloud base is generated by organized entrainment of air from levels below cloud base The fraction originating from each of these levels is set proportional to the level mass The lifted parcel is raised through the next model level During this process turbulent entrainment into Eu and detrainment out of D the cloud are parameterized
3. 21 21 21 23 23 24 24 24 25 25 25 25 25 25 25 25 28 30 31 31 33 35 36 41 42 45 45 46 47 47 17 Appendix 7 Flowtrace and other performance tests 53 17 1 TOMCAT Cray YMP8 2 ee 53 17 2 SUIMGAT lt Cray YMP sit oe EEE ag 1 AA a 54 reS SLIMGAT Cray YMPS s dd bo i a ia tated ft alla i eters Bat GS coo 55 ITA TOMCAT Green O3000 23 syr se a Re eRe ay eg dh BE NLR FR BE AR ey GG ew 56 1 Introduction This manual describes the TOMCAT SLIMCAT off line three dimensional 3D chemical transport model CTM The manual is aimed at providing both some background to the formulation approach of the model as well as some instruction on how to run it or at least understand the standard job decks Section 2 gives some background to the model Section 3 provides some basic definitions of atmospheric dynamics Section 4 describes the model grid Section 5 describes the model transport and advection Sections 6 and 7 describe the model physics and chemistry routines respectively Section 8 provides information on the sources of winds available to force the model Section 9 gives information on actually running the model Users simply interested in running the model could just focus on Sections 2 7 8 and 9 2 Background 2 1 History The TOMCAT and SLIMCAT off line CTMs have been widely used for atmospheric chemistry studies over the past decade or more The TOMCAT model was originally written in 1991 92
4. 1 20 Kise 2 7 rion B04 Rea EG ee tren 4 6 4 These functions Ka form a basis for the ensemble of the polynomial functions of order less than or equal to 2 over the box B In addition this basis is orthogonal because ies Ko x y Kg x y dedy 0 21 for all pairs a 8 in 0 z y z xy yy which satisfy a 4 6 However the basis is not normalised and we have X pY f I K3 a y drdy XY S 22 o Jo X pY f Ken ad 5 3 o Jo X pY K z y dedy 5 3 o Jo X pY K sa 5 5 o Jo X pY I Ky z y dzdy S 5 X Y 2 1 f Ky 7 9 dedy S 9 31 The moments of a function r x y over the box B are defined by the following equations Mo I j I REE a 23 X pY o Jo X pY M sf f p Ky x y r x y dedy o Jo X pY Maz 5 p Kg a y r 2 y dedy o Jo X pY Myy 5 p Kyy 2 y r a y dedy o Jo X Y Mszy of e p Koy 2 y x y drdy 0 0 where p represents the mass density taken to be uniform in the box B Equally knowledge of the moment of a function r x y allows us to have an analytic expression of the function or more exactly of its projection over the basis of Kg r x y X MaKa a y 24 where M pS is the total mass of the box B Interpretation of the moments As the functions K are dimensionless if the function r x y is the mass mixing ratio of a constituent the moments of r as defined above have the dimension of mass In pa
5. 1989 28 Williamson D L and P J Rasch Two dimensional semi Lagrangian transport with shape preserving interpolation Monthly Weather Review 117 102 129 1989 27 11 Appendix 1 Notation Symbol A B Ptop 3203 RH Su So SM k E SIIC AS S 3 Wim O Az Parameters to define a p levels radius of the earth Parameter to define g 0 levels heat capacity of dry air at constant pressure divergence Mass flux divergence Stability function Gravitational acceleration Constant for precipitation parameterization Asymptotic mixing length Vertical diffusion coefficient Mixing length Latent heat of condensation at 0 C Liquid water content Fraction of tracer mass from level j transferred through model interface k due to all sub gridscale processes Mass flux of air through interface k from updraft subsidence or vertical diffusion pressure surface pressure Pressure at top model interface Pressure at top of level Reference pressure Diabatic heating rate Specific humidity Gas constant Relative humidity Richardson s number Dry static energy Virtual static energy Tracer mass in level j Air mass in level j time temperature Surface temperature Mass flux of tracer through interface k from updraft subsidence or vertical diffusion zonal wind meridional wind u cos latitude v cos latitude horizontal wind vector Vertical velocity Vertical mass flux Verti
6. S Wofsy G Russell and D Rind Chemistry of the global troposphere fluorocarbons as tracers of air motion J Geophys Res 92 6579 6613 1987 26 19 Rood et al 1991 20 Russell G L and J A Lerner A new finite differencing scheme for the tracer transport equation J Appl Meteorol 20 1483 1498 1981 21 Shine K P The middle atmosphere in the absence of dynamical heat fluxes Q J R Meteorol Soc 113 603 633 1987 22 Shine K P Zonal momentum in the middle atmosphere Q J R Meteorol Soc 115 265 292 1989 23 Shine K P and J A Rickaby Solar radiative heating due to absorption by ozone in Proceedings of Quadrennial Ozone Symposium Gottingen pp 597 600 A Deepak Hampton Va 1989 24 Stockwell D Z and M P Chipperfield A tropospheric chemical transport model Development and vali dation of the model transport schemes Q J Roy Met Soc 125 1747 1783 1999 25 Stockwell D Z C Giannakopoulos P H Plantevin G D Carver M P Chipperfield K S Law J A Pyle D E Shallcross K Y Wang Modelling NOx from lightning and its impact on global chemical fields Atmospheric Environment 33 No 27 4477 4493 1999 26 Swinbank R and A O Neill A stratosphere troposphere data assimilation system Mon Weather Rev 122 686 702 1994 27 Tiedtke M A comprehensive mass flux scheme for cumulus parameterization in large scale models Mon Wea Review 117 1779 1800
7. SUBROUTINE FOR CHEMISTRY C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IMPLICIT NONE CALL PARADI CALL CSTUNI CALL FOR3D CALL REDEM CALL GRILLE CALL MOMENT C REAL FAC INTEGER I IAN IBIS JOUR JV K L MOIS NUMJOU C Work out day of year JOUR INT RJOUR MOIS INT RMOIS 18 IAN INT RAN CALL WCALEN NUMJOU IBIS JOUR MOIS IAN C Overwrite bottom level with current time DO 40 K 1 LAT DO I 1 LON SO I K NIV 1 RAN RMOIS 1 0 12 0 RJOUR 365 0 1 SM I K NIV ENDDO 40 CONTINUE Cc Cc Store advected tracers S0 in output arrays ST DO 50 L 1 NIV DO K 1 LAT DO I 1 LON ST I K L 1 S0 I K L 1 ENDDO ENDDO 50 CONTINUE RETURN END 7 5 GLOMAP GLOMAP is a set of subroutines which forms a module for calculating aerosols properties e g see Mann 2004 This can be included in the TOMCAT SLIMCAT model In this case the CHIMIE subroutine calls the GLOMAP aerosol subroutines as well as the ASAD chemical routines 7 6 DLAPSE DLAPSE is a lagrangian particle sedimentation scheme used to calculate denitrification in the polar strato sphere e g Davies 2004 Originally this code was a separate IDL model but it is now inserted into the TOMCAT SLIMCAT trajectory subroutines using F90 19 8 Forcing Winds and Temperatures This section describes the various sources of winds and temperatures that can be used to force TOMCAT SLIMCAT 8 1 ECMWF Analyses The ECMWF analyses
8. Thus for their experiments with a 3D off line CTM Prather et al 1987 used the following sequence for the advection steps e Advection along X At e Advection along Y At e Advection along X At e Advection along Z At e Advection along X At 34 e Advection along Y At e Advection along X At with At 4At and At 2Atz The validity of this approach depends on the shear of the wind rather than on its magnitude In the case of a uniform wind field the successive treatment of different directions does not introduce an error with respect to a simultaneous treatment however in regions of strong wind shear to treat directions successively can cause errors in the transport but it is difficult to estimate the effect quantitatively 13 3 Generalisation of the Scheme Case of a variable density In the above we assumed that the density was constant and uniform over all of the domain in fact the equations obtained remain valid even if this is not the case due to the property described above Suppose that we had calculated the moments of a box B from those of the two sub boxes with different density which compose it Considering box B of density pi and dimension X box B of density p2 and dimension X gt and box B of density p with _ Xipi X2p2 X We proceed as if the boxes B and B had a density of p therefore B becomes Bi of length Xi by compression and By becomes Bi of length X3 by stretchin
9. as Ey Mu u Dy My u where the fractional entrainment detrainment rates u u depend inversely on cloud radii Clouds of smaller sizes are assumed to occur in the absence of large scale convergence i e shallow convection and it is further assumed that 54 so that a simple form is used 1x107 m for deep convection Eu u 3x10 4 mt for shallow convection 2No test is made below cloud base for parcel buoyancy because it is assumed that large scale moisture convergence which is the final boundary condition for cloud occurrence sufficiently destabilizes the below cloud levels and ensures thermals reach the lifting condensation level LCL Note that if buoyancy is not achieved at the LCL in this parameterization scheme there is no further lifting of the parcel to test for a higher level of free convection 14 The updraft mass flux and tracer mass flux vary as OM og eae oTr Tr ae Eu Xe Du Xp Eu Xe Du Fp During ascent through the level T and q are adjusted along the wet adiabat and modified for entrainment of environmental air and detrainment of cloud air Any condensate is added to the liquid water loading Detrainment of liquid water to the environment is included At the level top the lifted parcel is tested for buoyancy as before if it is buoyant then lifting through subsequent levels continues otherwise the lower interface of the model level in which buoyancy is
10. at the Centre National de Recherches M t orologiques CNRM Toulouse and was first described and used by Chipperfield et al 1993 for studies of the polar stratosphere The TOMCAT model used hybrid o p levels Although it performed reasonably well it was not ideal for stratospheric studies and could not make best use of the stratospheric forcing analyses then available Therefore I developed the related SLIMCAT 3D CTM first described in Chipperfield et al 1996 This differed from TOMCAT in using an isentropic vertical coordinate and hence having a domain effectively limited to the stratosphere and using diagnosed heating rates to derive the vertical transport These developments allowed the model to use the UKMO analyses for multiannual stratospheric simulations e g Chipperfield 1999 Meanwhile the TOMCAT CTM was further developed for tropospheric chemistry studies by the inclusion for example of convection Stockwell and Chipperfield 1999 wet and dry deposition Giannakopoulos et al 1999 lightning Stockwell et al 1999 and detailed gas phase chemistry Law et al 1998 Effectively SLIMCAT became a stratospheric CTM and TOMCAT usually forced by ECMWF analyses which extended to only 10 hPa a tropospheric CTM Recent technical and scientific developments have meant that it would be desirable to combine the two related CTMs into a single model with a flexible vertical coordinate and different methods or treating key processes
11. condition imposed on the tracer distribution at the interfaces of the boxes e Upwind character in contrast to centred finite difference schemes advection by the Prather scheme of any feature only affects the grid boxes downstream of this feature CPU memory and time requirement The major inconvenience of the Prather scheme which explains why it is not more widely used is the cost in CPU memory and also CPU time For each advected tracer the scheme requires for N dimensions e 1 array for the zero order moment e N arrays for the first order moments e N N 1 2 arrays for the second order moments Thus for a 3D model using the conservation of second order moments there are 10 3D arrays required for each tracer and 6 2D arrays for a 2D model If the available memory is too limited for example in high resolution or with a large number of tracers it should be noted that it is possible to truncate the scheme to the conservation of first order moments only which limits the number of arrays per tracer to N 1 i e 4 in three dimensions Of course the accuracy of the scheme and the characteristic low diffusivity are affected by this operation according to Prather 1986 the scheme truncated in this way equivalent to the slope scheme of Russell and Lerner 1981 has an accuracy comparable to that of a 2nd order finite difference scheme However it should be noted that effective doubling of the resolution obtained by the
12. conservation of 2nd order moments can be achieved at a lower cost in CPU memory multiplication by a factor 2 5 in 3D and in CPU time depending on the problem than by increasing the number of gridboxes 13 4 Use in Spherical Geometry The Prather scheme first requires a grid to be defined over the model domain and second the calculation of the mass fluxes at the interfaces of all of the gridboxes thus defined Here we show how the scheme can be implemented in a transport model using spherical geometry and then how the mass fluxes at the box interfaces can be calculated exactly when the wind field is read in as spectral coefficients The results obtained here can obviously be used in an off line transport model forced by analyses or from output from a spectral model i e TOMCAT or alternatively for the on line transport of tracers e g water vapour in a spectral GCM Definition of the grid A rectangular grid in spherical geometry is defined by the latitudes of the north south interfaces and the longitudes of the east west interfaces i e by a series of latitudes fx k 0 1 K and longitudes Ai i 0 1 I The gridboxes Bir defined by Bir Ai 1 Ai x 0x 1 k cover the whole globe so that do 3 south pole ox north pole Ao r periodicity in longitude Note that at no stage do we assume that the grid is regular Any rectangular grid which satisfies the above conditions is suitable for use
13. correcting the first and second order moments before this advection to give a distribution which is positive everywhere the correction proposed by Prather 1986 consists of 3 3 Mz min 3 Mo maz Mo m 50 1 Mz min 216 gl Mz maz 1M Mo Mrz 51 Note that this correction slightly increases the diffusivity of the initial scheme but this effect is very small as the creation of negative values as described above is comparatively rare and the correction consists of a redistribution of tracer mass within a gridbox This correction is termed flux limiting as it consists of avoiding negative mass fluxes However the use of such a flux limiter will destroy tracer correlations Generally it is best to avoid the use of the limiter if negative tracer values can be tolerated 14 Appendix 4 Chemically Updating Tracer Moments When using the Prather 1986 advection scheme in the basic TOMCAT SLIMCAT model any change in tracer mass due to chemistry only acts upon the zero order tracer moment S0 Strictly the first and second order moments should also be modified This would increase the chemical resolution of the model and could reduce the rate of chemical reaction between two tracers with opposite gradients within a gridbox for example However chemically updating the lst and 2nd order moments would be prohibitively expensive in a full chemistry simulation Consider the integrals over a gridbox in two dimensions of th
14. e g vertical transport First in 1999 the ECMWF extended the upper boundary of their analyses to 0 1 hPa and recently completed the ERA 40 reanalyses from 1957 2002 over this same domain Second many key scientific issues concern processes in the upper troposphere lower stratosphere and so model boundaries in this region are not desirable This manual describes the new version of TOMCATi SLIMCAT which combine both original models and other routines into a single library The resulting model extends from the surface upwards with different options for the vertical coordinate The new library can still be used for all past applications of the earlier models with some significant new possibilities TOMCAT was originally an acronym for the Toulouse Off line Model of Chemistry and Transport though now it is just a name 2 2 Web Pages The home web page for the CTM is at www env leeds ac uk tomcat www env leeds ac uk slimcat Both addresses link to the same page On these pages are general information about the models an up to date list of references which have used the models with reprints of some papers links to results for recent campaigns etc There are also pages with restricted access for Users and Leeds based people Please get password from me 2 3 Nupdate and Libraries The TOMCAT SLIMCAT code is contained in a library which is managed using nupdate Nupdate is software originally developed for Cray supercomputers and pr
15. following equations the reverse of the above to calculate the moments of a box B from those of the two sub boxes BP and B which comprise it Mo My MP 28 M aMP 1 0 MP 3 1 a MP a ME Mar 0 MP 1 a ME 5a 1 a MP MS Maz 5 1 2a 1 a MP a MS My M MP My ME ME May aMP 1 a ME 3 1 a MP aMf where a M MP ME Reconstruction of a box When it is necessary to reconstruct a box from three sub boxes B Bz and B3 we proceed by first regrouping B and B using 19 then adding Bs to this result We can verify that the final result is identical to that obtained if we begin by first regrouping Bj and B3 and then add B The same result is valid for the splitting of a box into 3 sub boxes Alternate Directions Time splitting In a multi dimensional motion the method of evaluating the total advection by successively calculating the advection along each direction is known as time splitting This method is not the only one possible because in principle it is possible to proceed in one step by splitting each box into 3N sub boxes according to the destination of the tracer mass However it is is much easier and less costly to proceed by separating the directions In addition this method gives the possibility of using a different timestep according to the direction in relation to the different Courant numbers in each direction
16. lost is set as the top of the cloud column ktop see Figure 2 Centre of Levels Level Interfaces cr Ktop 2 Keio Du gt Dyer EE EEE NN ktop 1 ktop Du gt Dyr laap ktop Figure 2 Showing cloud top and detrainment of cloud air Organized detrainment occurs above cloud top and allows for some small shallow cumuli overshooting the level of zero buoyancy Tiedtke 1989 It is parameterized as Die 1 8 Mir Diep 1 B M ter 1 ge 0 3 for shallow convection 0 0 for deep convection Discretization PME ATi At a model level k ME ME g tt p The tracer updraft mass flux can be implicitly calculated as Tr u ME Trt Pett Bit yt DEN 3No modification is made to the large scale environmental thermodynamic properties qe Te even though a change in these would be expected due to detrainment of cloud air 4Note that there is turbulent entrainment into and detrainment out of level ktop in addition to organized detrainment 15 Substituting M MEt ER DEF and So t SM 1 and solving for Tr gives sott DEE Tru Trat But gan y gt u u 19 The updraft tracer flux through interface k can be represented as a linear function 1 Trt V 500 op j N where j is the fractional tracer mass from level j transferred through interface k in the convective updraft From the condition that at the lowest model interface M MN Tr 0 it c
17. uniformly spreading the mass flux over the whole length of the interface is equivalent to substituting the real wind field with a discretised one which has the following properties e The mass flux across the interfaces of the grid are the same with the discretised wind field as with the real wind field In particular if the real wind field is non divergent the same is true of the discretised wind field e The meridional velocity v is uniform over the north south interfaces e The zonal velocity u cos is uniform over the east west interfaces In addition it follows that the discretised wind field has discontinuities at the interfaces of the gridboxes there is a discontinuity in the meridional component at the east west interfaces and a discontinuity in the zonal component at the north south interfaces The problem of the treatment of the pole a singularity in a gridpoint model is a classical problem A test to verify the performance of a numerical advection scheme in the case of cross polar flow is to consider a 40 solid rotation about and axis inclined at 90 degrees to the polar axis This test has been used for example by Williamson and Rasch 1989 for a semi lagrangian scheme taking an initial tracer distribution of a localised structure at the equator of the solid rotation The equation of the motion is the following u UsinAcosg v UsinxX As the motion is a solid rotation the initial structure moves on the spher
18. 455 Total samples 494 550 Accumulated time secs 10 0 Time per sample msecs 4 Sample bin width bytes index secs cum samples function dso file line 1 112 710 22 8 22 8 11271 pris a out prog f 44844 2 69 030 14 0 36 7 6903 harwel a out prog f 26656 3 45 620 9 2 46 0 4562 __mp_barrier_nthreads libmp so mp_barrier c 271 4 37 240 7 5 53 5 3724 ftoy a out prog f 35354 5 27 780 5 6 59 1 2778 jac a out prog f 42894 6 23 950 4 8 64 0 2395 impact a out prog f 38164 7 23 650 4 8 68 7 2365 __exp libm so exp c 103 8 13 700 2 8 71 5 1370 __mpdo_advx2_1 a out prog f 9071 9 12 620 2 64 74 1 1262 diffun a out prog f 34321 10 12 210 2 5 76 5 1221 __mpdo_consom_1 a out prog f 12777 11 11 970 2 4 79 0 1197 __pow libm so pow c 256 12 11 860 2 4 81 4 1186 __mpdo_advy2_1 a out prog f 9683 13 10 350 2 1 83 4 1035 acaceto a out prog f 28479 14 9 410 1 9 85 4 941 __mp_barrier_master libmp so mp_barrier c 363 15 8 070 1 6 87 0 807 __mpdo_advz2_1 a out prog f 10228 56 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 lt snip gt Oo 00000000 00000r RRP KK KK NDNNPDNDPDND000 00 494 850 840 190 130 060 900 530 470 370 250 980 610 490 330 190 050 010 910 8
19. 5296 qpassm a out prog f 15875 presetr a out prog f 78043 pharwel a out prog f 26080 __mpdo_wdepin_4 a out prog f 21008 __Mmp_wait_for_completion libmp so mp c 883 ludcmp2 a out prog f 17760 calsub a out prog f 12858 quanto13 a out prog f 29156 vdiff a out prog f 84934 __mpdo_physics_2 a out prog f 70943 minv a out prog f 17673 TOTAL 57
20. 7 CLOUD 1 94E 00 32768 5 91E 05 1 10 94 87 FINCYCL 1 80E 00 4 4 51E 01 1 02 95 89 QSAT 1 38E 00 292132 4 71E 06 0 78 96 67 DQSATDT 8 62E 01 168175 5 13E 06 0 49 97 16 LOUIS 6 44E 01 32768 1 96E 05 0 36 97 53 PELF 5 85E 01 5 1 17E 01 0 33 98 22 SUBSCAL 5 05E 01 32768 1 54E 05 0 29 98 51 REEMDT 4 55E 01 1 4 55E 01 0 26 98 76 PEFL 4 51E 01 160 2 82E 03 0 26 99 02 PEFP 3 27E 01 320 1 02E 03 0 19 99 20 TOMCAT 2 63E 01 2 63E 01 0 15 99 35 INIEXP 2 58E 01 2 58E 01 0 15 99 50 REEZNOT 2 49E 01 2 49E 01 0 14 99 64 PEPF 2 46E 01 320 7 70E 04 0 14 99 78 ADVEC 1 67E 01 48 3 49E 03 0 09 99 88 INICYCL 7 08E 02 4 1 77E 02 0 04 99 92 CEP 3 25E 02 320 1 01E 04 0 02 99 93 CEF1 3 04E 02 320 9 51E 05 0 02 99 95 CONVEC 2 33E 02 48 4 85E 04 0 01 99 96 ALRET 2 27E 02 5 4 54E 03 0 01 99 98 53 CEF2 1 26E 02 320 3 92E 05 0 01 99 98 CORPOLE 8 40E 03 5 1 68E 03 0 00 99 99 CHTRON 7 03E 03 385 1 83E 05 0 00 99 99 FINITER 6 79E 03 48 1 41E 04 0 00 100 00 INICSF 2 64E 03 5 5 27E 04 0 00 100 00 WRCHK 2 12E 03 1 2 12E 03 0 00 100 00 FINEXP 8 28E 05 1 8 28E 05 0 00 100 00 INICSTE 3 46E 06 1 3 46E 06 0 00 100 00 INCHK 2 26E 06 1 2 26E 06 0 00 100 00 Totals 1 76E 02 725361 Jun 12 19 31 mpc 77 95 Mflops 191 93s a out The most expensive routines are the second order moments advection scheme and CONSOM which applies the convection diffusion matrix to the tracers However the code is written efficiently for the Cray 17 2 SLIMCAT Cray YMP8 Below is a flowtrace from
21. 70 810 750 730 720 680 660 660 620 600 460 440 420 550 100 o00O0O0O0O0O0O0O0O0o00o00000000000000000000Nxe 2 2 6 6 6 6 5 5 5 5 4 3 3 3 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 0 88 89 90 90 91 91 92 92 93 93 94 94 94 95 95 95 95 95 96 3 96 96 96 96 97 97 2 97 97 97 97 96 97 100 2 3 0 6 2 8 3 8 3 8 2 5 8 1 3 5 Th 9 1 4 6 Th 8 0 1 3 4 5 6 0 585 584 319 313 306 290 253 247 237 225 198 161 149 133 119 105 101 91 87 81 75 73 72 68 66 66 62 60 46 44 42 49455 radabs a out prog f 72194 __powdi libftn so pow_di c 49 fyinit a out prog f 36795 asad_step a out prog f 32437 __read libc so 1 malloc c 907 compiled in read s memcpy libc so 1 stat c 32 compiled in bcopy s __mpdo_chimie_1 a out prog f 18622 _defgu2sd libffio so defgu2s c 243 lubksb2 a out prog f 17728 __mpdo_calflu_17 a out prog f 4876 log libm so log c 207 Write libc so 1 flush c 58 compiled in write s fyfixr a out prog f 36286 fft991 a out prog f 15642 _get_value libfortran so lread c 735 calflu a out prog f 3925 rpassm a out prog f 16625 spetrul a out prog f 79519 _ld_read libfortran so lread c 239 calcpv a out prog f 1
22. ASAD chemical solver UNIASADx x listing in uniasadx x list and a library for the CCM boundary layer scheme UNIPBLx x listing in unipblx x_list 2 4 Fortran The original TOMCAT and SLIMCAT models were written in quite standard Fortran77 On the basis that there is no point rewriting something which works fine these routines are still Fortran77 Some newer additions to the model library use Fortran90 e g the chemical data assimilation code the DLAPSE microphysical model Depending on the components selected for a particular run then a f77 compiler may be ok An f90 compiler will be ok in all cases Some notes on writing Fortran code for TOMCAT SLIMCAT are given in Appendix 6 2 5 Parallelisation TOMCAT SLIMCAT was originally written to run efficiently on vector machines e g old Cray machines This is achieved by making inner DO loops vectorisable and the code should still conform to that However the high performance computers HPCs available for running the model now are parallel machines There are different methods for parallelising code These methods can be divided into distributed memory or shared memory categories In distributed memory machines a copy of the code is run on each processor and when necessary messages are passed between the processors to communicate data This approach requires that the code be specifically written for the parallel machine using for example Message Passing Interface MPI In shared me
23. ATMX are similar lengths 15 3 Compiling and Running As for normal model get an example job deck 15 4 MPI Performance The execution time T of a parallel program defined as the time elapsed between the first processor beginning and the last processor completing the task depends on a number of variables The principal ones are the size of the problem N grid resolution and number of layers and tracers in this case and the number of processors P Execution time is often not the most effective way of assessing performance Since T varies with problem size even on a single processor it is difficult to compare performance at different N Efficiency E removes this size dependence and is defined as E T PTp 47 Table 1 CPU time speedup and efficiency for a single tracer at T42 NPE time s speedup efficiency 1 fuji 9483 1 T3E 7530 1 00 1 00 2 3700 2 04 1 02 4 1900 3 96 0 99 8 943 7 99 1 00 16 504 14 98 0 94 32 296 25 44 0 80 64 213 34 35 0 54 128 192 39 23 0 31 Table 2 CPU time speedup and efficiency for MPI job with a single tracer at T42 NPE time s speedup efficiency initialisation pbl transport chemistry total compared to 4PE fuji 4704 4 65 270 1395 3595 5325 8 49 138 1242 1853 3281 6 49 0 81 16 23 77 741 875 1716 12 41 0 78 32 7 65 491 393 955 22 30 0 70 64 9 43 425 231 708 30 09 0 47 where Tp is the execution time on P processors Speedup S is the factor by which execution time is reduced on P p
24. LT scheme 24 9 6 Output Files The model produces a range of output files and brief details are given here Usually the user ftps the output files from the remote machine which runs TOMCAT SLIMCAT to his her local machine for analysis and storage Most of the output files are single precision 32 bit unformatted binary files Reading them on some machines i e linux will require software that switch from big endian to little endian binary format T have a selection of fortran and IDL programs to read these files I also have jobs which convert the output e g PDG files to NetCDF files 9 6 1 PDG The PDG file is the basic output from the model run It contains the global 3D fields of tracers temperature PV humidity and pressure at a requested frequency controlled by NSO1 9 6 2 DIAG This file contains basic information on the model transport terms e g mass fluxes winds etc Mainly used for debugging diagnosing the model Controlled by NS02 9 6 3 STAT This file contains profiles of model tracers temperature pressure and PV interpolated to a specified set of stations currently set to 55 This allows frequent sampling of the model at just a few sites Controlled by NSO3 9 6 4 ZON This file contains zonal mean output on the 15th day of each month This is a smaller filer than the PDG file and gives a more regular sampling for multiannual runs 9 6 5 OXO 9 6 6 REST This file is for restarting a new TOMCAT SLIMCAT run
25. RT TRUE flowview Luc Flowtrace Statistics Report Showing Routines Sorted by CPU Time Descending CPU Times are Shown in Seconds Routine Name Tot Time Calls Avg Time Percentage Accum ADVX2 2 73E 01 480 5 70E 02 22 92 22 92 or ADVY2 2 28E 01 480 4 76E 02 19 15 42 07 CALFLU 1 08E 01 21 1 13E 00 19 00 61 07 INCIRA 1 40E 01 21 6 68E 01 11 76 73 83 ADVZ2 1 32E 01 480 2 74E 02 11 03 84 69 HRTCO2 5 03E 00 1344 3 75E 03 4 22 88 91 SOLAR 2 46E 00 1344 1 83E 03 2 06 90 97 HRTO3 2 17E 00 1344 1 61E 03 1 82 92 79 INITER 2 07E 00 240 8 61E 03 1 73 94 52 FINCYCL 1 41E 00 20 7 04E 02 1 18 95 71 INTRAN 1 30E 00 1 1 30E 00 1 09 96 80 HTRH20 8 10E 01 1344 6 03E 04 0 68 97 47 ADVEC 7 84E 01 240 3 27E 03 0 66 98 13 SETUSS 7 61E 01 1344 5 66E 04 0 64 98 77 INIEXP 2 85E 01 1 2 85E 01 0 24 99 01 FINITER 2 72E 01 240 1 13E 03 0 23 99 24 CHIMIE 2 38E 01 240 9 92E 04 0 20 99 44 INICYCL 2 17E 01 20 1 08E 02 0 18 99 62 THERML 1 11E 01 1344 8 27E 05 0 09 99 71 PCMDRD 9 76E 02 21 4 65E 03 0 08 99 79 CHTRON 7 54E 02 4071 1 85E 05 0 06 99 86 SOLANG 6 68E 02 1344 4 97E 05 0 06 99 91 SLIMCAT 6 03E 02 1 6 03E 02 0 05 99 96 PCMRAD 3 84E 02 1344 2 86E 05 0 03 99 99 CORPOLE 5 65E 03 21 2 69E 04 0 00 100 00 ORBIT 2 25E 04 21 1 07E 05 0 00 100 00 WCALEN 8 23E 05 21 3 92E 06 0 00 100 00 CALNUM 8 18E 05 1 8 18E 05 0 00 100 00 FINEXP 8 15E 05 1 8 15E 05 0 00 100 00 GENGRID 3 55E 05 1 3 55E 05 0 00 100 00 INICSTE 3 64E 06 1 3 64E 06 0 00 100 00 Tota
26. The Semi Lagrangian Scheme The SLT scheme implemented in TOMCAT SLIMCAT is based on the method of Williamson and Rasch 1989 and closely follows the details given in their paper 10 5 3 1 Extended Grids The SLT scheme requires the wind velocities and tracer mixing ratios to be known on a grid which extends two gridpoints beyond the standard model grid in the X Y and Z directions This is required by the cubic interpolation scheme Cyclic continuity is used in longitude In latitude the grid is extended to include a pole point row and one row across the pole In the vertical linear interpolation is used in the upper and lower two levels to avoid the need to create points outside the model domain 5 3 2 Interpolation The interpolation method has an important impact on the performance of a semi Lagrangian advection scheme In the scheme described here the wind velocities and tracer mixing ratios are interpolated in each of the X Y and Z directions using a cubic Lagrange scheme 5 3 3 Coordinate Transformation For trajectory positions polewards of DLATLIM see below the SLT scheme performs a coordinate transformation to a local geodisic coordinate The local geodisic coordinate is essentially a rotated spherical coordinate system whose equator goes through the original poles This is based on the transformation described by Williamson and Rasch 1989 except only one transformation is performed rather than a different transformation for e
27. The TOMCAT SLIMCAT Off Line 3D CTM User s Manual Martyn Chipperfield Institute for Atmospheric Science School of Earth and Environment University of Leeds Volume I A Description of the Model and a Guide to Using it Draft Still being written August 2006 Version 0 8 Contents 1 2 Introduction Background 21 History dea SIA E TT ED aly ek el E AE Me 22 Web Pages ii ean AA Sr Sr re Ge on ne Co ee ee ee Gen Rte eo 2 32 Nupd te and Libraries y uke a Gite tated at AU Nede hel kat Se eee 24 Fortrano Gece ele Sl gel sh ref dessa ka eet ie ea hb ft Judy eA 2 0 ParallelisatiGm ec ra SUS A kle ee Se eae A PE 2 5 1 OpenMP Parallelisation ee 2 5 2 MPI Parallelisation rv vr vr vr vr rn rn E E E 2 6 Platforms so 23 e022 a a sce EE Ee ee bedrer nr Tee es ka adel d Basic Dynamics Model Grid 4 1 Horizontal grid lt a DA ee eR BAA OR Ne ar ar SR ee la 42 Vertical Grid x va ge eh er BPS ey A A Ee ee eae 65 Transport and Dynamics orl Mass FlUxES asc Bl IS Ad Sok Sri oh we ge et Sheth let SEE eee Ee a 5 1 1 Consistency of Vertical Horizontal Mass Fluxes 2 aar vr vr rv vr rv ev 5 1 2 Consistency of Advected Mass and Analyses 2 vr vr vr e 5 2 Advection The Prather Scheme 2 2 0 ee ee 5 3 Advection The Semi Lagrangian Scheme 2 2 rv rv rv ee ee Bra Extended Grids ii A dd Loud DUS SEE eee we ee De 5 3 2 Interpolation 2 shes es do de d
28. a 10 day single level run with NIV 1 LON 128 LAT 64 with 24 hourly UKMO winds and LVERT FALSE flowview Luc Flowtrace Statistics Report Showing Routines Sorted by CPU Time Descending CPU Times are Shown in Seconds Routine Name Tot Time Calls Avg Time Percentage Accum ADVX2 1 86E 01 960 1 94E 02 32 77 32 77 Add kkk ADVY2 1 53E 01 960 1 60E 02 26 95 59 72 eR CALUVI 4 66E 00 11 4 24E 01 8 20 79 53 CHTRUV 2 60E 00 11 2 36E 01 4 57 84 10 CHTRTG 2 26E 00 11 2 06E 01 3 98 88 09 CONUKM 1 71E 00 11 1 56E 01 3 02 91 10 INITER 1 53E 00 480 3 19E 03 2 70 93 80 ADVEC 1 51E 00 480 3 15E 03 2 66 96 46 INICYCL 1 20E 00 10 1 20E 01 2 11 98 58 FINCYCL 2 35E 01 10 2 35E 02 0 41 98 99 SLIMCAT 1 56E 01 1 1 56E 01 0 28 99 27 FINITER 1 53E 01 480 3 19E 04 0 27 99 53 INIEXP 1 38E 01 1 1 38E 01 0 24 99 78 CHIMIE 7 49E 02 480 1 56E 04 0 13 99 91 CALFLUK 5 17E 02 11 4 70E 03 0 09 100 00 FINEXP 8 27E 05 1 8 27E 05 0 00 100 00 CALNUM 8 05E 05 1 8 05E 05 0 00 100 00 GENGRID 3 49E 05 1 3 49E 05 0 00 100 00 INICSTE 3 35E 06 1 3 35E 06 0 00 100 00 Totals 5 08E 01 6321 The 2 Prather advection routines x and y directions are the most expensive The next 3 most expensive routines CALUVI CHTRUV CHTRTG simply interpolate the UKMO winds to the SLIMCAT grid These could 54 probably be made more efficient 17 3 SLIMCAT Cray YMP8 Below is a flowtrace from a 5 day multilevel run with NIV 3 LON 128 LAT 64 with 6 hourly ECMWF winds and LVE
29. ach trajectory point 5 4 Other Advections Schemes WAF 5 5 Trajectory Calculations TOMCAT SLIMCAT can calculate Lagrangian air parcel trajectories The model can either be used for just this purpose or Langrangian trajectories can be included in a normal Eulerian TOMCAT SLIMCAT run for comparison The same options for trajectories exist as for the main CTM e g vertical coordinate method for vertical motion The trajectory code uses either a simple Euler forward time scheme or by default a more sophisticated 4th order Runge Kutta time scheme The Runge Kutta scheme is based on that of Fisher et al 1993 The switch ITRAJ in the subroutine EVOL determines which timestepping scheme is used 11 6 Physics 6 1 Convection The convection scheme implemented in TOMCAT SLIMCAT is based on the work of Tiedtke 1989 It is a mass flux scheme which explicitly represents the cloud fields and their circulation in order to determine the effect of convection on the large scale budgets of heat and moisture Representing this circulation is vital for tracer transport because the convective induced fluxes are needed Other schemes such as the convective adjustment scheme of Betts and Miller 1984 and the Kuo 1974 scheme do not explicitly represent this circulation By using a simple bulk model to represent the cloud ensemble the contribution of convection to the large scale budget equations of heat and moisture can be calculated in the Tiedt
30. an be seen that By 0 j 1 2 N i e there is no flux through the lowest model interface From equation 19 a recursion formula can be used to determine general column values of PL k 1 k 1 E DE Peg Pra Unit ado mM gett where 5 E de if j k 1 adjacent levels k j 1 O otherwise 6 1 2 Subsidence Outside of the cloud a sub gridscale subsidence flux is induced to maintain mass balance within the column The mass flux through interface k MF is given by MI Mi is Trk Mk yk E so that the tracer flux is Tr M Xe and OL Smr 04 3 16 6 2 PBL Louis 1979 or Holtslag and Boville 1983 6 3 Radiation Within SLIMCAT a radiation scheme is used to diagnose heating rates to calculate the vertical motion in the isentropic part of the model There are different options for the radiation scheme Temperature Domain Longwave Shortwave Albedo Clouds Refs Details of CTM Radiation Schemes MIDRAD Analyses 0 700 hPa CO 15 um O3 9 6 um vibro bands H20 Oz and O3 6 bands 125 175 175 205 nm 206 244 244 278 nm 278 363 408 853 nm 0 2 Ay gt 278 nm None Shine 1987 Shine and Rickaby 1989 17 CCMRAD Analyses 0 hPa surface O3 CO H20 CHy N20 F11 F12 0 3000 cm AA 100 cm Oz 03 CO H20 18 intervals 200 5000 nm ZHMRAD Analyses O hPa surface Climatology fn Aw t O Aw Not included Briegleb 1992a b 7 C
31. ant t Considering a box B figure 1 and motion along the direction Ox 3 cases are possible depending on the values of U and U the mass fluxes across the interfaces Bi 1 Bi and Bi Bi 1 a U and U are of the same sign b U positive and U negative convergence c U negative and U positive divergence The box B is this split into a maximum of 3 sub boxes BF contains the air which is going to leave the box B by the left face if U is negative BP contains the air which is going to leave the box B by the right face if Uj is positive B contains the air which is going to stay in the box B during the timestep 2 Calculate the sub moments of each sub box The following equations allow the calculation of the sub moments of two sub boxes BP and B of the box B split into two along Oz For the right sub box BP X X x 0 Y MP a Mo 1 a Mz 1 a 1 2a Mzz 26 MP M 3 1 a Mzz MP o Moe M Mzy MP a M 1 a May M aMyy and for the left sub box BY 0 X x 0 Y ME 1 a Mo aM a 1 2a Mzz 27 MS 1 a Mz 3a Mzz 33 M 1 a Ma ME 1 a May ME 1 a My aMzy Mg 1 a Myy where a X X X represents the fraction of mass of the box B contained in the box BP 3 Regroup the different sub boxes to reconstruct the new boxes at the instant t dt We use the
32. ations between x 4 U V and y C are the following 1 OU OV PO NESSE ER 1 4 zm lat S Sa 1 OV OU Se EEE i a 1 p Es l a n Vx 32 C Vw _ Ox 2y OW Ow 2 Ox In the following we will also use the zonal 7 and meridional n divergences defined by n n n 34 37 A 1 ow T Faw Or 35 1 ov Pir 2 ees VE 5 36 It is possible to calculate exactly the mass fluxes at the interfaces of a rectangular grid on the condition that this grid is regular in longitude but not necessarily in latitude of which a Gaussian grid is a particular example 1 Zonal mass flux Let FA be the mass flux across the interface between the boxes B and B y1 p By definition Pr Fi f puad 37 k 1 which can also be written ia i F f p gt Ud Hk 1 1 u 2 Meridional mass flux Let F be the mass flux across the interface between the boxes B and Bi k 1 By definition Ai pia f pvacos d 38 Ai 1 which can be written de Ff 1 pVdr A 4 1 If we consider the box B our problem is therefore to calculate FA p F s Fig and F y a calculation of the meridional flux The first step consists of calculating the Fourier coefficients of V at the north south interfaces of the grid by using the values of Pam tabulated at the latitudes of the interfaces For example Vm 10 Vi Yi 39 n Performing a direct Fourier transform on these coefficients will give us the values of V at
33. cal mass flux at top interface Geopotential height Vertical depth of the model level 28 6 371x10 m 1005 46 J K kg kg m 257 adimensional ms 2 0 x 107 438 18 m m s7 m J kg kg water kg air kg air m7 s Pa Pa 100000 Pa K s kg water kg air 287 05 J kg adimensional J kg J kg kg tracer kg air K kg tracer m7 s7 Pa s kg m s kg s Ratio R Cp Tracer mass mixing ratio latitude longitude 3 Von Karman s constant Wavelength sing Fraction of tracer mass from level j transferred through interface k from updraft subsidence or vertical diffusion UA as S Density kg m7 Water density kg water m Tracer density kg tracer m Air density kg air m7 relative vorticity vertical coordinate Potential temperature Reference 0 pressure coordinate vertical velocity Superscripts and Subscripts Model level k 29 12 Appendix 2 Universal Constants The following table lists the variable names of the universal constants in TOMCAT SLIMCAT They are contained in the fortran common block CSTES which is in the nupdate common deck CSTUNI Variable Dewiption 0 CONV CP CPSG CPSL CPV CPVMCP DEUOMG ECPH ETV ETVQ GG GSCP GSRA OMEGA RA RASCP RASCP2 RASL RASRV RTER RTER2 RV RVSRA STEFAN TMERGL TOO UNSCP UNSG VKARMN XL XLF XLI XLISCP XLISG XLSCP XLSG XLSRV XPI XPOO 180 XPI Specific heat capacity of dry air at co
34. ction with nupdate common decks or included files so that code is only in 1 place e Don t insert nupdate common decks in a SUBROUTINE unnecessarily 16 3 Example Subroutine The following subroutine illustrates this style SUBROUTINE EXAMPLE AIN BIN COUT NDIM CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c C Calculate lt gt c C Method Cc c lt description gt C c Input C C AIN lt description units gt C BIN lt description units gt C c Output Cc C COUT lt description units gt C C References C c lt details of any papers reports etc gt C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 50 IMPLICIT NONE c CALL PARADI CALL MOMENT CALL FOR3D c C Arguments INTEGER NDIM REAL AIN NDIM BIN NDIM COUT NDIM c C Local arrays REAL SUM INTEGER I K L N c c c Calculate mass weighted mean global T K SUM 0 0 DO L 1 NIV DO K 1 LAT DO I 1 LON SUT SUT T3D I K L SM I K L SUM SUM SM I K L ENDDO ENDDO ENDDO SUT SUT SUM FLOAT LON LAT NIV c C Scale DO N 1 NDIM COUT N AIN N BIN N SUM ENDDO c RETURN END 16 4 Example Interface A model generally consists of different modules i e collection of subroutines to calculate different processes in the atmosphere e g chemistry radiation microphysics Although these modules need to be coupled the code in each module should be largely indepen
35. ctory model UGAMP internal report no 37 November 1995 9 Fisher M A O Neill R Sutton Rapid descent of Mesospheric air into the Stratospheric Polar Vortex Geophys Res Lett 20 1267 1270 1993 10 Giannakopoulos C M P Chipperfield K S Law and J A Pyle Validation and intercomparison of wet and dry deposition schemes using Pb 210 in a global three dimensional off line chemical transport model J Geophys Res 104 23 761 23 784 1999 11 Heimann M The global atmospheric tracer model TM2 Technical Report no 10 Deutsches Klimarechen zentrum Hamburg 1995 12 Holtslag A A M and B Boville Local versus nonlocal boundary layer diffusion in a global climate model J Clim 6 1825 1842 1993 13 Kuo H L Further studies of the parameterization of the influence of cumulus convection of large scale flow J Atmos Sci 31 1232 1240 1974 14 Law K S P H Plantevin D E Shallcross H L Rogers J A Pyle C Grouhel V Thouret and A Marenco Evaluation of modeled O3 using Measurements of Ozone by Airbus In Service Aircraft MOZAIC data J Geophys Res 103 25 721 25 737 1998 15 Louis J F A parametric model of vertical eddy fluxes in the atmosphere Boundary Layer Meteorol 17 187 202 1979 16 Mann GLOMAP 2004 17 Prather M J Numerical advection by conservation of second order moments J Geophys Res 91 6671 6681 1986 18 Prather M J M McElroy
36. dent Indeed it is an advantage to have a module mostly self contained as it will be more portable and can be more easily coupled with other models In TOMCAT SLIMCAT many parts of the model are completely entwined for example the code to read the analyses calculate mass fluxes output the results This is really the core of the model and makes large use of nupdate common decks to pass variables around However when adding new modules it is better if they interface with TOMCAT SLIMCAT in just 1 routine The module can still use nupdate common blocks to store pass its own variables but it is good to avoid using TOMCAT SLIMCAT ones Variables can be passed as arguments If a few constants are needed rather than calling the TOMCAT SLIMCAT deck CSTUNI the constant can be declared as a local PARAMETER and given the same name value SUBROUTINE INTERF ol CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c C Example interface routine C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IMPLICIT NONE c CALL PARADI CALL MOMENT CALL FOR3D c INTEGER NDIM NAER PARAMETER NDIM LON NIV NAER 2 C C Aerosol fields REAL T1D NDIM P1D NDIM S01D NDIM S01D NDIM NAER c Cc Read in fields on LAT LON and NIV grid here c Cc Loop over latitude can be parallelised DO 100 K 1 LAT c c Copy TOMCAT 3D arrays for e g tracers T p to C arrays for aerosol calculation DO L 1 NIV DO I 1 LON II L 1 NIV LON T
37. e oe eek Bee eRe tte a RE te aa KE 5 3 3 Coordinate Transformation 5 4 Other Advections Schemes ee 5 5 Trajectory Calculations a Physics 6 1 Convection Saar ERROR D Gud ale geet eens ale SNE TETTE tah SA de NS ar 6 1 1 Convective Updraft ee 6 1 2 Subsidence a ee eae a ee eae ea eS DR 62 PBL a a eee ee ee ih bes fre EGG dressene 485 63 Radiation a eras aly PET A GE ee ee ee eat thane da AAA 4 Chemistry Schemes TE VOtratOSPherie or ve ces Bethy BOR AYO a are se Ede EET Ee kG 7 2 ASAD based Tropospheric ee des Canolle OZone son ag A a eed ia oe ween RS NA EEE ft e E TA User Supplied EA ee E ITA A A AAA ee Ta GLOMAP A a A iaa RE Dd a a TE DEABSES A A eta He wy A aaa Ss 10 10 10 10 10 10 11 11 11 11 11 12 12 12 16 17 17 8 Forcing Winds and Temperatures Sal EE MNP Analyses n E PEE EA EE Oww E S EORR A SESS SEES ewe 8 2 UK Met Office UARS Analys g scoici kr a e Pw ee ae Oe Bo EO 9 Running the TOMCAT SLIMCAT Model OT P rameters sad Alt Oe be dex Six fag fad ara a O tte Ain OS a A Ke E a tiara En 9 2 Model Switches and Variables ee 9 3 Fortran Channels zd orar Do Sabel A A A SSS SG gAs Subroutines Lac a ai E Gu al Gi Glace wg daa dd OD NAVECEION A a a o a AE S i Seb A ee NE dte Bh nee OU a RR 0 5 1 Prather chee a Feat OE ee EE e 05 2 SLI Scheme is AA A A A EEE FG AAA es 96 Output Biles got tele dt A ADAN bd oe 961 PDGA waa a o aan e
38. e products of any 3 basis functions of the moments equation 11 above If one of the functions is the zero order moment K we get X Y f i Kok3 x y drdy XY S 52 0 0 X Y i KoK z y dedy 3 53 0 0 A Y 2 P Kok x y dedy S 3 54 xX Y f KoK 2 y dedy 8 5 55 0 0 X Y 1 KoK3 x y dedy S 5 56 0 0 X Y f Koki 2 9 azdy 5 9 57 0 0 42 Otherwise the integrals give Ko y dr 0 58 rato en de 0 59 REG a Ky a ae 25 45 60 T K x y de 25 35 61 Koen k ed ken dedy 5 9 62 a K2 0 y Kes 2 0 dedy 28 45 63 Similar results are obtained for K and Kyy or K and K All of the other products involving the cross terms mixtures of x y and z are zero see equation 12 above Consider two chemical tracers whose distributions in two dimensions are given by A x y 5 Y AaKa x y 64 and TB z y SV Bak z y 65 The rate of the chemical reaction between the two species depends on the product of the concentrations rate k 4 B rate k ApKo AsK AzoKgo AyKy AyyKyy A0y KoyMBoKo B K2 BooKoo ByKy ByyKyy BoyKoy The effect on the zero order moment is proportional to Er 1 1 1 de ABKodxdy AoBo 3 A B AB 5 Agr Bez AyyByy gav Bey 66 So in addition to the product of the average concentration of with the average concentration of B first term on RHS the zero order moments of and B are modifi
39. e within the vertical column is maintained by including sub grid subsidence of environmental air induced by convection within the same timestep Subsidence is parameterized differently to Tiedtke 1989 6 1 1 Convective Updraft The N vertical levels in TOMCAT extend from level M closest to the Earth s surface to level 1 e g at 10 hPa in the stratosphere with old ECMWF grid Thermodynamic and other variables are defined either at the centre of model levels or at level interfaces see Figure 1 Air at the top of level M i e at interface M 1 is lifted through model level NV 1 along the dry adiabat qT 9 dz Cp The notation used in this report is summarised in Appendix 2 The temperature of the lifted parcel at the top of level M 1 i e interface M 2 is given by PA o dn gjen 12 Centre of Levels Updraft and Level Interfaces Subsidence Variables O skred orde 0 1 s Pa SEN 1 2 PS sn te 2 vane Oe mites RAMs cues ati Rt ans k 1 k EFS Dt gt ee MEA ME oo k kt FS ES EN MEA MEH kven k l k 2 ed eh je k 2 N 2 date FE NR VE gg ND N 1 EE TNL REN ak MT N ter N Figure 1 Convections variables shown on TOMCAT SLIMCAT vertical levels Although the specific humidity remains constant with parcel lifting at this higher level the temperature of the lifted parcel is lower than at interface NV 1 and therefore the relative humidity RH increases If the lifted parcel becomes supersa
40. e without being deformed and returns to its initial position after travelling around the globe and crossing the two poles A method sometimes employed to solve this deficiency is to cover the pole by a polar box a circular zone centred on the pole inside which the mixing ratio of tracers is taken to be uniform Therefore the treatment of the singular point is avoided This method has been used by Rood et al 1991 with the Van Leer scheme and seems capable of giving good results Therefore it was tested with the Prather scheme by taking the polar box to be the ensemble of boxes which touch the pole However this approach gave transport across the pole which was too rapid Therefore in TOMCAT the mass flux across the pole is obtained by calculating the wind vector at the pole and determining the mass of each triangular box which would be transported to the diametrically opposite box in each timestep This then gives the mass flux across the pole 13 5 Correction of Negative Values As advection schemes are not perfect it is necessary to couple them with a certain number of fixes to correct their intrinsic deficiencies A well known problem with classical transport schemes spectral and finite difference is their non positiveness For minor constituents e g water vapour ozone it is obvious that negative values are purely a numerical artefact and that these negative values can be a problem in chemical or physical calculations Contrary to what
41. ed by products of the higher order moments The effect on the first order moment in the x direction is proportional to a 1 2 1 f ABK dady 5 Az Bo AoBz 15 Ac Bre Are Br 9 AyBay AzyBy 67 o Jo and similarly for the y direction The effect on the second order moment in the x direction is proportional to X Y 1 2 2 2 J ABK dxdy zz Bo AoBaz Az Bz Azz Baz AzyBay 68 o Jo 5 15 35 45 43 The effect on the second order moment in the xy direction is proportional to X Y ABK dxdy AoBry AzyBo Ay Bz Az By Z AzyBaz Aro Boy AyyBay AzyByy 69 Therefore to chemically update the 0 1st and 2nd order moments of a chemical tracer in three dimensions would require integrating the chemistry 10 times As chemical integration is generally by far the most costly part of a CTM this is prohibitively expensive However it may be practical to use the higher order moments to update just the zero order moments equation 57 44 15 Appendix 5 MPI version of TOMCAT SLIMCAT As described the Section 2 5 the current default version of the model is parallelised using OpenMP This will run on many platforms although it is possible that in the future hardware issues may limit high end HPC machines to distributed memmory parallelisation Therefore I have created a copy of the standard library UNICATMPPO0 xx which is identical to the standard library as far as possible but wi
42. f mass passing directly from one polar gridbox to the gridbox diametrically opposite the pole is therefore not always zero The mass flux from one box to the other occurs by crossing one by one the neighbouring boxes therefore it is a case of a zonal mass flux and not a meridional one Unfortunately the zonal flow to a neighbouring box is poorly represented near the pole in the Prather scheme as discussed above In effect one of the characteristics of the scheme is that the mass flux calculated at the interfaces of the boxes are spread over these interfaces as of they were due to a uniform velocity this does not mean that the calculation of the mass flux does not take account of the wind shear on the interface Thus in planar geometry the volumes transported from one box to another are rectangular slabs the velocity being assumed uniform over the sides of the boxes in our case of interest spherical geometry the volumes transported along longitude are in reality angular sectors which is equivalent to considering the angular velocity a dd dt u cos constant over the interface In particular for the boxes neighbouring the pole the motion is that of a rotation about the pole which returns to taking a zero velocity at that point In summary the scheme implicitly considers in part that the mass flux across the pole is zero which is true but equally that the velocity at the pole is zero which is not true in general In fact the operation of
43. f processors to be used PARAMETER NPROCMAX 64 PARADI 76 NSLSTMX 3000 PARADI 77 NRLSTMX 3000 PARADI 78 NFLDSL 220 PARADI 79 NSLSTMX specifies the maximum number of grid points any processor will need to send information about to any other processor in preparation for the advection step and NRLSTMX is the equivalent buffer size for receiving data The total number of model fields to be exchanged prior to advection is given by NFLDSL which is set to MAX 3 NIV 1 10 NTRA for Prather SOM advection For most choices of processor number there are several possible grid decompositions and the possibilities can be calculated using the program calc_slcons f which is in the home martyn UNICAT directory The parameters in the program need to be sent to match the actual TOMCAT SLIMCAT job for LON LAT NIV NTRA NDYN DTO and UMAX The program then ouputs the possible domain decompositions e g 46 Partial output from calc slcons f Possible choices for SLCONS are Number of processors 64 nlonmx 8 nlatmx nhalmx nproci nprock 8 nprocmax 64 nslstmx 48 nristmx 48 nfldsl 8401 nlist nlonmx 16 nlatmx 2 nhalmx 8 nproci 4 nprock 16 nprocmax 64 nslstmx 56 nristmx 56 nfldsl 8401 nlist 2 Number of processors 32 The model generally runs quicker if the patches are as square as possible as this minimises communication between processors so choose a decomposition where NLONMX and NL
44. from the final time of the previous run 10 Acknowledgements I am grateful to Pascal Simon of the CNRM Toulouse for many helpful discussions Others 25 References 1 Andrews D G J R Holton and C B Leovy Middle Atmosphere Dynamics Academic Press 1987 2 Bacmeister J T M R Schoeberl M E Summers J R Rosenfield and X Zhu Descent of long lived trace gases in the winter polar vortex J Geophys Res 100 11 669 11 684 1995 3 Betts A K and M J Miller A new convective adjustment scheme ECMWF Technical Report No 43 1984 cs Briegleb B P Delta Eddington Approximation for Solar Radiation in the NCAR Community Climate Model J Geophys Res 97 7603 7612 1992 5 Carver G D P D Brown and O Wild The ASAD atmospheric chemistry integration package and chemical reaction database Computer Physics Communications 105 197 1997 D Chipperfield M P D Cariolle P Simon R Ramaroson and D J Lary A three dimensional modelling study of trace species in the Arctic lower stratosphere during winter 1989 90 J Geophys Res 98 7199 7218 1993 7 Chipperfield M P M L Santee L Froidevaux G L Manney W G Read J W Waters A E Roche and J M Russell Analysis of UARS data in the southern polar vortex in September 1992 using a chemical transport model J Geophys Res 101 18 861 18 881 1996 8 Chipperfield M P J Kettleborough and A Pardaens The TOPCAT offline traje
45. g if pa gt pi this does not change their moments Thus it is the case of regrouping two boxes of identical density and this is achieved by using 19 Extension to three dimensions The extension of the definitions and of the equations 17 18 and 19 to three dimensions is straightfor ward for example the equations 17 must be completed by D 2 MP Ms D MP aM D MP aM Properties of the scheme The Prather scheme possesses the following properties e Stability the scheme forward Eulerian is stable in the limit a lt 1 which is rigorously equivalent to the CFL condition udt dx lt 1 e Conservation the formulation of the scheme in flux form assures the exact conservation of the total mass of each tracer e Accuracy compared to finite difference schemes for which the mixing ratio is only known at the centre of the gridbox additional information is provided by the storage of the first and second order moments and the calculation of their evolution Prather 1986 estimated that the conservation of second order moments conferred on his scheme an accuracy comparable to a fourth order finite difference scheme 35 e Small diffusivity and local character the scheme is well adapted to the representation of localised phe nomena because the advection is performed by exchange between adjacent boxes therefore there are no distance effects with this scheme It should be noted that there is no continuity
46. g the following nupdate commands I deck xx Insert after line deck xx B deck xx Insert before line deck xx D deck xx yy Delete the range of lines deck xx to deck yy and insert any lines which follow The user can also control which decks are included in the model Fortran file This avoids the need to include and compile every line in the library in every job and is also a way in which alternative parameterisations can be selected The nupdate decks are selected by using C Therefore part of a typical TOMCAT SLIMCAT job may start cat lt lt eof gt tomi up C SLIMCAT INIEXP INICSTE INICYCL INITER FINITER FINCYCL FINEXP C ALRET CALFLU CORPOLE CHTRUK CHTRON CALNUM RADIAT lt lt select other decks needed gt gt IDENT EXP D PARADI 10 11 C Model grid PARAMETER LON 64 LAT 32 NIV 31 C Number of tracres PARAMETER NTRA 30 NVTOT 48 lt lt other temporary updates gt gt eof The official version of Cray nupdate does alot more than this but it is also very difficult to install on a system Therefore a stripped down version written in plain f77 specifically for the CTM is available This has a limited set of the nupdate commands and is easy to use on any platform it just needs a Fortran compiler and makes the model easily portable Overall there are 3 nupdate libraries relevant for TOMCAT runs In addition to the TOMCAT SLIMCAT library itself there is a library containing the
47. hemistry Schemes A range of different chemical schemes can be included in TOMCAT SLIMCAT These can be standard full chemistry schemes for the stratosphere or troposphere or can be simple schemes written supplied by the user In all cases the basic approach is that the TOMCAT SLIMCAT job has a chemistry subroutine called CHIMIE If this subroutine is short and does not modify the model tracers then they are passive On the other hand the CHIMIE subroutine can call the subroutines for a full chemistry scheme which may include emissions deposition etc as well as the integration of the chemistry itself 7 1 Stratospheric There is a standard stratospheric chemistry scheme contained in the TOMCAT SLIMCAT library in the nupdate deck TOMCHI See separate manual Volume II for details 7 2 ASAD based Tropospheric The full tropospheric chemistry scheme is based on the ASAD solver Carver et al 1997 which is contained in its own nupdate library UNIASADx x The standard code for the default scheme is contained in the nupdate deck TROPCHI 7 3 Cariolle Ozone parameterised O3 7 4 User Supplied The user can supply his her own version of the CHIMIE subroutine to modify the advected tracers to suit any particular experiment An example of a subroutine for a single age tracer is given here SUBROUTINE CHIMIE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C c DUMMY INTERFACE
48. iD II T3D I K L P1D II PL I K L SO1D II 1 50 1 K L 10 S01D 11 2 S0 I K L 11 ENDDO ENDDO c C Call aerosol scheme CALL AEROSO T1D P1D S01D NDIM NAER C C Copy aerosol results to TOMCAT tracers DO L 1 NIV DO I 1 LON II L 1 NIV LON SO I K L 10 S01D II 1 SO I K L 11 S01D II 2 ENDDO ENDDO c End of parallel latitude loop 100 CONTINUE c RETURN END 52 17 Appendix 7 Flowtrace and other performance tests Below is a flowtrace from a 1 day multilevel SLIMCAT run using the CCM radiation scheme with NIV 13 LON 64 LAT 32 with 24 hourly UKMO winds and LVERT TRUE This can be compared with an equivalent run using MIDRAD 17 1 TOMCAT Cray YMP8 Below is a flowtrace on a Cray YMP8 from a 1 day run with NIV 19 LON 128 LAT 64 with 6 hourly ECWMF winds and 2 tracers The job included convection and vertical diffusion flowview Luc Flowtrace Statistics Report Showing Routines Sorted by CPU Time Descending CPU Times are Shown in Seconds Routine Name Tot Time Calls Avg Time Percentage Accum ADVX2 4 03E 01 96 4 20E 01 22 88 22 88 kkkk ADVY2 3 40E 01 96 3 54E 01 19 28 42 15 ero CONSOM 2 73E 01 48 5 69E 01 15 49 57 64 ADVZ2 1 59E 01 48 3 32E 01 9 02 66 67 CALSUB 1 31E 01 4 3 28E 00 7 43 74 10 CONVMA 9 88E 00 32768 3 01E 04 5 60 79 70 MUHERM 5 76E 00 131072 4 40E 05 3 27 87 58 CHIMIE 4 74E 00 48 9 87E 02 2 69 90 26 CALFLU 3 80E 00 5 7 59E 01 2 15 92 42 INITER 2 39E 00 48 4 98E 02 1 36 93 7
49. ing resolution cp MARTYN UTIL TRONxx fort 2 files at model grid resolution cp MARTYN UTIL TRONyy fort 1 cp MARTYN UTIL LEGCyy fort 20 cp MARTYN UTIL LEGIyy fort 21 The above data files need to be correctly set according to the values of MIO and MI1 xx MII and yy MI0 The model outputs a number files On channel 25 the model writes a PDG file which contains the 3D tracer arrays at intervals specified by the variable NS01 On channel 31 the model outputs REST file which can be used to restart a continuation run jour n results mv fort 25 EXP UNIOO1 PDGO1 mv fort 12 EXP UNIOO1 DIAG1 mv fort 15 EXP UNIOO1 ZONO1 mv fort 31 EXP UNI001 REST1 mv fort 14 EXP UNIOO1 TRAJ1 4 The following lines in the subroutine INIEXP control the length of the model timestep and frequency of output DTO 3600 0 INIEXP 56 NITERT 24 INIEXP 57 NDYN 1 INIEXP 58 NS01 24 INIEXP 59 NS02 24 INIEXP 60 NS03 24 INIEXP 61 NS04 24 INIEXP 62 NS05 1 INIEXP 63 DTO is the basic model timestep This is split into NDYN dynamical subtimesteps NITERT is the number of iterations in one cycle i e the time between the forcing files divided by DTO Output is written to the PDG file every NSO1 iterations to the DIAG file every NSO2 iterations to the STAT file every NS03 iterations and to the TRAJ file every NSO4 iterations NSO5 is a switch to write to the ZON files 1 or not 0 The Prather advection scheme can use a limiter to prevent negative mi
50. ion of the levels follows equation 11 For the hybrid levels the definition of the model levels changes with altitude Above a reference potential temperature o the model uses pure isentropic levels and the potential temperature of an interface is defined as 0k43 Cho Between o e g typically 350 K and the surface the model uses hybrid levels The pressure of the model half levels in this region is then given by Prot C poo Dps where poo is the pressure at the lowest purely isentropic half level and D 1 C This method of defining a vertical coordinate is possible in an off line CTM as p at the surface ps and at o are known at future times using ps and T from the analyses Bottom boundary 5 Transport and Dynamics A key part of any off line CTM is the routines which convert the wind temperature analyses used to force the model to the CTM grid and the tracer advection schemes The conversion of the wind temperature fields should be done in such a way that the analyses are not degraded or modified unduly The tracer advection method should ideally have the desirable properties of conservation low numerical diffusion etc This section first describes the way in which TOMCAT SLIMCAT maps the analysed wind fields onto the model grid The following subsections described the main tracer advection routines available in the model 5 1 Mass Fluxes 5 1 1 Consistency of Vertical Horizontal Mass Fluxes 5 1 2 Co
51. is sometimes assumed the formulation of the Prather scheme does not ensure the posi tiveness of a tracer without an additional correction This can be seen in the example in figure 3 of Prather 1986 In a one dimensional example if the tracer at to is localised in a single box B of the grid with the 1st and 2nd order moments set to zero and the wind u is directed towards the right M 1 46 Mi 0 Mi 0 Mt 0 47 MH 0 Mit 0 TT After a timestep dt the structure is globally displaced towards the right and only the distributions of the boxes B and B 1 are modified the distribution at ty dt in the box B 41 calculated using 17 18 and 19 is the following Mj 1 a 48 Mi 3a 1 a Mi 5a 2a 1 1 a Mi a 49 Mit 3a 1 a 41 Mit 5a 2a 1 1 0 At the following timestep for certain values of the wind a the quantity of tracer mass which will be transported into the box Bi 2 will be negative which bearing in mind the absence of tracer in this box at to dt will create a negative values at the instant to 2dt As can be seen the appearance of negative values is due to the fact that during the reconstruction of the boxes at to dt using the formulae 19 a distribution is created which is globally positive Mo gt 0 even in the box B 1 but not the local positiveness allowing the advection of negative quantities of mass Therefore this problem can be solved by
52. ke scheme as a v Vs og a Fi Musa Masa Mu Ma 3 L cu ea amp p Sr AC Qr 17 ty V t ge Mugu Maga Mu Majd Ueu ea amp amp FS Gun 18 where s is the dry static energy q the specific humidity Pair the density of air v the horizontal velocity vector w the vertical velocity c the rate of condensation e the rate of evaporation and Qpr the radiative heating Overbars denote grid area means primes denote deviations from the mean is the evaporation of cloud air that has been detrained into the environment and is the evaporation of precipitation in the unsaturated subcloud layer The index tu denotes boundary layer turbulence For both expressions the first two terms on the right side represent the perturbation to the heat and moisture budgets by convection The convection scheme implemented in TOMCAT SLIMCAT is identical to that described by Tiedtke 1989 apart from the following differences in TOMCAT midlevel convection and convective downdrafts are not included and there is no organized entrainment of environmental air above cloud base The scheme does include cumulus updrafts in the vertical column entrainment of environmental air into the cloud and detrainment of cloud air to the environment The magnitudes of these are related to horizontal convergence of moisture below cloud and the difference between cloud and environmental specific humidity at cloud base Mass balanc
53. ls 1 19E 02 17417 Again the Prather advection schemes account for about 53 of the total cost The subroutine INCIRA is relatively expensive as it just interpolates the CIRA Oz and T fields onto the SLIMCAT model levels In general the MIDRAD scheme is well coded and cheap Overall the above job ran at 72 Mflops 55 Considering the number of calculations involved in the Prather second order moment advection scheme the subroutines ADVX2 etc are coded efficiently These are the most computationally intensive routines and should dominate the total cost of the model As can be seen from the flowtraces above they do The rest of the model is also sufficiently efficient Of course when a chemistry scheme is coupled to the model this will then dominate the total cost 17 4 TOMCAT Green 03000 Below is a profile from a 1 day TOMCAT full tropospheric chemistry run on the CSAR 03000 Green 16 CPU The run had NIV 31 LON 64 LAT 32 with 6 hourly ECMWF winds prof lines a out pcsampx m SpeedShop profile listing generated Mon Jul 26 11 35 50 2004 prof lines a out pcsampx m a out n64 Target program pcsampx Experiment name pc 4 10000 0 cu Marching orders R12000 R12010 CPU FPU 500 Number of CPUs 400 Clock frequency MHz Experiment notes From file a out pcsampx m Caliper point 0 at target begin PID 10765979 tmp jtmp tmpdir_126594 a out Caliper point 1 at exit 0 Summary of statistical PC sampling data pcsampx 49
54. m outside This is done by creating a halo of grid points around each patch which the processor is aware of but does not perform any computations on This is illustrated in figure 4 for processor 12 of 32 on a grid of 8 LAT x 16 LON In order to keep the data associated with these halo points up to date the processors need to communicate with their neighbours before each advection time step The size of halo required depends on the length of the time step and the resolution of the grid since it must contain all points from which tracers may enter the domain during the time step The time taken to perform the communication depends both on the size of the message to be communicated and on a startup time during which the processors set up a communication channel Thus communication time increases as the time step increases and more importantly as the number of processors increases The latter occurs for two reasons Firstly more processors means more halos to be filled and secondly communication may be required not just between adjacent processors but also between more remote neighbours This can result in large overheads and a consequent decrease in the efficiency of the program Communication is achieved in UNICATMPP via the message passing interface MPI 45 In principle it is possible to parallelise the code to allow for irregular non rectangular patches of different sizes Cate Bridgeman s orginal code allowed for this in certain but not all ro
55. mory machines the code is run on all the processors which work on a subset of the common memory This approach is usually much easier to implement in existing codes e g using OpenMP 2 5 1 OpenMP Parallelisation The standard version of TOMCAT SLIMCAT has been parallelised using OpenMP In OpenMP work is gener ally shared at the DO LOOP level As such existing codes can be parallelised in stages e g so that only the most expensive loops are shared over all the processors and the rest of the code is just run on a single processor the Master thread The disadvantage of OpenMP is that the efficiency of the parallelisation of individual loops may not scale that well with increasing numbers of processors To get good parallel speed up the parallelisation should be done at the highest level possible e g the chemistry is parallelised over an outer latitude loop which contains all the routines to solve the chemistry on a longitude xlevel slice With this the parallelisation speedup for the chemistry an expensive part of the model is effectively 100 efficient The speedup is similar for other expensive parts of the model such as the radiation schemes and boundary layer scheme which have an outer latitude loop which can be easily parallelised In effect it is only the dynamical subroutines where the speed up is not as efficient A consequence of this method of parallelisation is that the number of CPUs used to the run the model should be a facto
56. nsistency of Advected Mass and Analyses In SLIMCAT the total mass is an advected quantity After 1 advection step this mass should ideally be equal to that diagnosed from the analyses used to force the model However as the vertical mass fluxes are diagnosed from a radiation scheme and due to spatial temporal interpolation errors this will not be the case It is possible to correct the mass fluxes within SLIMCAT so that Mass above a theta level psup l plt 1 ggdphidlambda 12 The mean mass flux over the cycle top wgrs l 0 5 wgri wgriz dphidlambda 13 Mass change above model interface over cycle dm psup psup 14 Correction to integrated mass flux corr dm dt 15 which is spread evenly over surface dwgri corr 4arr 16 5 2 Advection The Prather Scheme The default advection scheme used in TOMCAT SLIMCAT is a gridpoint scheme which is based on the con servation of second order moments of the tracer distribution within the gridboxes This scheme possesses the properties of being stable accurate conservative and is one of the best performing of the transport schemes available The second order moments scheme was developed by Prather 1986 The model can also be used with the conservation of first order moments equivalent to the slopes scheme of Russell and Lerner 1981 or conservation of only zero order moments The Prather advection scheme is described in detail in Appendix 3 5 3 Advection
57. nstant p CP GG CP XL Specific heat capacity of water vapour at constant p CPV CP 2 QMEGA CPV CP 1 RVSRA 1 1 RASRV Acceleration due to gravity GG CP GG RA Earth s speed of rotation gas constant for dry air RA CP RA 2CP RA XL RA RV radius of the earth RTER RTER gas constant for water vapour RV RA ice melting temperature 1 CP 1 GG von Karman constant Latent heat of condensation at 0 C XLI XL latent heat of sublimation XLI CP XLI GG XL CP XL GG XL RV T reference pressure 30 Value 1005 46 JK kg 1869 46 JK kg 9 80665 ms 7 292x109 rad s7 287 05 JK kg 6371229 m 461 51 JK kg 5 6697x10 8 ms 271 23K 273 16K 0 4 2 5008x106 Jkg 2 83456x10 Jkg 3 14159 10 Pa 13 Appendix 3 Details of the Prather Advection Scheme After describing the principle of the advection method in the simple case of 2D transport in cartesian coordinates we will show how it is used in the more general 3D case in spherical coordinates 13 1 Definitions Consider a domain D in the cartesian grid Oxy covered by a regular grid of spacing X and Y in the following discussion the boxes of the grid will be referenced by the indices n in the direction Ox and k in the direction Oy Considering a portion B of the grid of which for simplicity the lower left hand corner is the origin of the grid i e B 0 X x 0 Y we can define the local functions on this portion Ko x y
58. ovides a nice environment for maintaining a centralised standard version of a model while still letting individual users make temporary changes in their own experiments Personally I think nupdate is still the best and easiest way of keeping control over large programs With nupdate the basic model library exists as a series of decks and within each deck every line of model code has a unique line number When you run a job you combine selected decks from the standard library with temporary changes to individual lines and nupdate produces the final Fortran program which can then be compiled The current TOMCAT SLIMCAT library is contained in the file UNICAT0 8 This file cannot be edited and so there is also a listing file unicat0 8 1ist which can be viewed and used to identify line numbers see martyn BIBLI unicat0 8_1ist An example of the listing file is COMDECK PARADI PARADT 1 C PARADI 2 INTEGER ITRIGO ITRIG1 IIFAXO IIFAX1 LAT PARADI 3 INTEGER LATS2 LEV LIGMAXO LKGMAXO LKG2 LOLA LON LPRO PARADI 4 INTEGER MIO MI1 MLON MLAT MPO MP1 PARADI 5 INTEGER NIV NMPO NMP1 NMIO NMI1 NPAR NPVAR NTRA NVRJPJ NVTOT PARADI 6 C PARADI 7 C parametres principaux du modele PARADI 8 C PARADI 9 PARAMETER LON 64 LAT 32 NIV 2 PARADI 10 PARAMETER NTRA 25 NVTOT 35 PARADT 11 PARAMETER NVRJPJ 0 PARADI 12 PARAMETER LOLA LON LAT LPRO LON LAT NIV LATS2 LAT 2 PARADI 13 The user can make temporary changes to the standard library by usin
59. quivalent to the CFL condition uA Az lt 1 for advection along Ox While this constraint does not pose any particular problems and allows the use of reasonable timesteps for the advection along latitudes or in the vertical it is not the same for the advection along longitude because of the convergence of the meridians towards the poles and the reduction of Ax which follows Ax acos AA To avoid very small timesteps the model uses extended polar zones as named by Prather et al 1987 at high latitudes This correction consists of grouping several adjacent boxes situated on the same latitude circle and advecting this block along the longitude in the same way as a single box of larger size The number of boxes grouped together depends on the latitude on the size of the boxes AA on the timestep At as well as on the maximum wind that will be encountered The number of boxes grouped together at each latitude is NUM which is set up in the routine CALNUM Before the advection step in the x direction the boxes are regrouped by using 17 then when the zonal advection is complete the original boxes are reconstructed using the reciprocal formulae 19 Problem of the flux across the pole As shown above the meridional mass flux Fg across an element of surface situated at latitude is given by Ai Fe f pva cos d dA 45 i 1 In particular the mass flux across the pole is zero even if the cross polar wind is not zero The quantity o
60. r of the number of model latitudes LAT for optimum efficiency For a scalar non OpenMP job the OpenMP directives inside the code are ignored 2 5 2 MPI Parallelisation Appendix 5 describes an MPI version of TOMCAT SLIMCAT which is available for use on HPC machines which do not support OpenMP or for which MPI jobs can use more processors e g on HPCx OpenMP jobs are limited to 16 CPUs The MPP library is equivalent to the standard library and allows the major options to be run e g full stratospheric or tropospheric chemistry runs 2 6 Platforms In the past TOMCAT SLIMCAT has been successfully run on a wide range of machines e g Cray 2 YMP J90 Hitachi Fujitsu Sun workstations SGI workstations and linux PC including my laptop Currently the default machines for the model runs are the SGI 03000 machines at CSAR Wren Green Fermat At HPCx the libraries and standard jobdecks for TOMCAT SLIMCAT are located in hpcx home n02 n02 1rkd unicat There are readme files which list the files there 3 Basic Dynamics This section lists some basic equations of atmospheric motion which are used in the detailed description of the transport part of TOMCAT SLIMCAT Section 5 The symbols used are summarised in Appendix 1 The continuity equation is Op oe _ Op Op On Ot Bay t aq 57 0 1 V vh oi Oy Integrating this equation from the top of the atmosphere to the ground gives the rate of change of surface pressure
61. rameters in TOMCAT number of latitudes number of longitudes spectral truncation of winds after conversion spectral truncation of winds read in number of latitudes in gridpoint forcing file number of longitudes in gridpoint forcing file number of levels 9 2 Model Switches and Variables The following switches are contained in the common deck SWITCH and can be set in the model jobdeck Meaning ICAT Vertical grid 1 SLIMCAT sigma theta eee 2 TOMCAT sigma pressure Vertical motion 0 None 1 Divergence 2 Mixed Divergence heating rates 3 Heating rates ICONV Convection 0 none 1 Tiedtke scheme Vertical mixing 0 none 1 Louis scheme 2 Simple complete mixing 3 Holtslag and Boville Dry deposition 0 none 1 Old scheme 2 New scheme needs IVDIF 3 Wet deposition 0 none 1 Giannakopolous Scheme 1 2 Giannakopolous Scheme 2 3 Giannakopolous Scheme 3 The length of the TOMCAT model run is controlled by the variables read in on channel 94 cat lt lt eof gt fort 94 0 NO OF FILES TO JUMP IN FORCING FILE 0 0 INITIALLY 1 RESTART 40 NCYCLT NDAYS NO FORCING FILES DAY 0 NFFILE NO FORCING OUTPUTS PER FILE eof The first line and last line can be ignored The second line indicates whether the start data comes from an initial data file 0 or from a restart file 1 The third line specifies the length of the run in numbers of analysis times 21 file at forc
62. roach to that for the calculation of the meridional mass flux The meridional divergence is obtained immediately by the difference between the meridional fluxes at the north and south interfaces of the box Finally we obtain the divergence of the zonal mass flux which allows us to calculate in turn the zonal flux itself at the east west interfaces of the gridboxes in terms of a constant from the condition of periodicity in longitude This constant is determined from the zonal mean of the zonal mass flux This method requires the prior tabulation of the values of the Legendre functions averaged over the boxes in the employed grid also in the file LEGIyy Note that it would be possible to calculate directly the zonal flux without using the divergence by using the formula A ON OE 42 1 p2 1 2 0A Op but this requires the additional tabulation of the derivatives of Pym This allows the calculation in terms of a constant from the periodic condition over the longitudes This constant is determined from D F An D F An 43 i i i e the calculation of the zonal mean of the zonal wind from the values at the centres of the boxes or from the values at the interfaces leads to the same result Thus the total mass flux F entering a box Bj is given by Fix Foi Ft Pees E and we can verify that Fik Sik Nik 44 High latitudes 39 The stability of the Prather scheme is conditional on the constraint lt 1 rigorously e
63. rocessors S PEp Table 1 shows CPU time in seconds speedup and efficiency for a 30 day run at T42 and 31 levels using a single methane tracer simple chemistry and Louis vertical diffusion scheme For this run speedup increases up to 64 processors but by this point communication costs outweigh the advantage of each processor dealing with fewer grid points and speedup is not substantially greater than that achieved using half as many processors Efficiency is nearly 100 for up to 16 processors and then begins to fall as input output overheads and communication costs become more important Table 2 shows CPU time in seconds for a 2 day run at T21 and 31 levels with NTRA 27 and NVTOT 49 using the PBL scheme and ASAD chemistry scheme 32 processors is a sensible choice for a run of this sort In moving from 32 to 64 processors the speedup increases by 35 but efficiency drops to below 50 It is apparent from the above times that the bottleneck occurs in the transport step There are two reasons for this Firstly it is the transport step that requires most of the communication between processors Secondly and more importantly quantities from the ECMWF forcing files undergo Fourier transforms before they are used in TOMCAT These transforms must be performed on an entire latitude and so a T21 run with 32 latitudes can only share the work between 32 processors Any additional processors will be idle while waiting for a transform to be completed and thi
64. rticular the zero order moment Mo represents the total mass of tracer contained in the box B The moment of order 1 in a given direction represents the average gradient over the box in this direction the first order moments also give the position relative to the centre of the box B of the centre of mass of the tracer whose distribution is represented by r x y The second order moments are proportional to the curve of the tracer distribution and give the matrix of the moments of inertia of the tracer distribution r x y for the box considered The variance of the distribution inside the box can equally be determined from the different moments TL pr x y dxdy Mo 25 X pY 1 1 1 f I pr x y dady Mo 3 Ma My 5 Mos Myy May 0 0 An important property which arises from the above definitions is that the moments of a function r z y contain the intrinsic properties of the grid in the sense that they remain invariant under all transformation on the coordinate system in other words all solid translations and or stretches compressions of the box B do not change the moments 32 Bi B Bi 1 Figure 3 Example one dimensional grid 13 2 Principle of scheme Using the above definitions the advection of a tracer over the period t to t dt is performed in several steps 1 Split each box of the domain into several sub boxes as a function of the mass flux at the interfaces between the boxes at the inst
65. s results in loss of efficiency 48 16 Appendix 6 Writing Fortran for TOMCAT SLIMCAT TOMCAT SLIMCAT are written in Fortran and in a certain style By following a few rules and conventions any new code should be easier to understand tidier more efficient and with a reduced risk of containing errors This document explains these rules and conventions Section 1 describes good practice for all Fortran programs Section 2 gives information specific to TOMCAT SLIMCAT 16 1 General Rules in Fortran 16 1 1 When Writing Code When writing new Fortran code you should always Declare all variables explicitly Use IMPLICIT NONE at the start of each subroutine However it is still good practice to use variables beginning with I N for INTEGERs and other letters for REALs Declare all constants as PARAMETERs Avoid using variables of 1 character except for loop indices Simplify mathematical expressions as much as possible before coding them up Avoid very short subroutines Reduce the use of intrinsic functions such as EXP power X Y etc which can be expensive by tidying up expressions and storing and reusing part expressions Read input files only once and store numbers if necessary Avoid large levels of nested loops 16 1 2 When Testing Code Get rid of all compiler warning messages Compile the code using different compilers to help get rid of non standard features and make it portable Use the compiler option to se
66. sd eek eo I AD RR FOR GOE Pe do f fe a 962 2DIAGK da E ns hh a a co at tM eae 963 STAT A A esker NT 964 LON ee SE See eee a EE are 065 CORO ok eek eg terse we ae ae hae ke ste Kok PS dar are BE TIERS 2 96 67 REST v re bier Se ae ee SE GE RE ah 10 Acknowledgements 11 Appendix 1 Notation 12 Appendix 2 Universal Constants 13 Appendix 3 Details of the Prather Advection Scheme VSL Definitions Avd at Be ee a G Re HORS g Aus EEE mogo doge ADN 13 21 Principle Of schemes sus 25 18 min ko hode hg OP stellet det dt AR A O Asti oh ok oe el 13 3 Generalisation of the Scheme 2 2 ee 13 4 Use in Spherical Geometry 13 5 Correction of Negative Values ee 14 Appendix 4 Chemically Updating Tracer Moments 15 Appendix 5 MPI version of TOMCAT SLIMCAT 15 1 Domain Decomposition 0 00 vr renn ee 15 27 Defining the Grid vas sad ee a oP Harada AR BR eS 15 3 Compiling nd Running sok kok eek A ee ee a 15 4 MPI Performance 16 Appendix 6 Writing Fortran for TOMCAT SLIMCAT 16 1 General Rules in Fortran eve vr arr ee ee 16 11 When Writing Code xx sv sh oe vg og ag S de Ga ll At RUE RR A OS OO ES 16 1 2 When Testing Code 4 are ST s r ba ene Hee dre ak SD le Ed 16 2 Styles Conventions in TOMCAT SLIMCAT o ooo rv rv vr vr rv nevn 16 3 Example Subroutines desea tee by AE dd Ao EGG GE fa a 16 4 Example Interface sui ltd pt HE RLE Le Sn AGRONOM leie elle GA 20 20 20
67. t all variables as undefined e g frt Hu Test the array bounds and consistency of subroutine calls e g frt Haes However this usually has a large overhead so turn this off for real experiments The most important point is to write correct code Once this code is running there is software e g ssrun on Green Wren Fermat which shows where the model is spending its time Then the code can be speeded up and optimised taking care not to affect the results 49 16 2 Styles Conventions in TOMCAT SLIMCAT Most of the code in TOMCAT SLIMCAT follows a certain style Given the size of the existing code compared to anything you may add it will be better if your new code fits in with this style Exceptions to the generic style occur if whole modules of code are coupled with the model For example the CCM radiation scheme and ASAD is written in a different style with lower case variables and so where these variables are referenced within TOMCAT lower case is also used I would say the conventions that define the TOMCAT SLIMCAT style are e The code is written in upper case The comments are written in normal text i e mixture of upper lower case e Code within loops and in IF ENDIF blocks etc should be indented 2 spaces e Three dimensional arrays are declared as LON LAT NIV The standard loop indices for these 3 dimen sions are I K and L respectively JV is used for tracers e Only use COMMON blocks in conjun
68. th MPI parallelisation This MPI version is based on older versions of parallel TOMCAT and SLIMCAT developed by Cate Bridgeman This appendix describes things specific to this MPI library Not all of the UNICAT options can be used under MPI The most likely use of the MPI library is to perform otherwise time consuming runs Therefore the MPI library is largely restricted to the default options e g Prather second order moments advection and the most expensive configurations of the model e g full chemistry 15 1 Domain Decomposition EE ojoje s 1617 101020202 for 2 20 20 fi 163775 rs 29 aoa aet 20350 oo rs 22300 aa e 24 25 25 20 20 27 27 2628 20 20 50J30 aa Figure 4 Computational grid for MPI model showing the domain owned by processor 12 and a possible halo 12 3 4 5 6 ae 10 11 12 13 14 15 16 Pppp 0000 DEE s ste fo fo jojos 4 10 8 The MPI library has been parallelised in the Single Program Multiple Data SPMD style which means that the same program is run by each processor but using its own data This is achieved in UNICATMPP by splitting the horizontal grid into patches so that each processor owns a unique subset of grid points and performs computations just on its own patch Sometimes a processor will need to know about points owned by neighbouring processors This is most important for the advection step when the processor needs to work out how much of a particular tracer has moved into its patch fro
69. that can be used to force TOMCAT are stored on HPCx in the directories hpcx home n02 n02 1rkd data ECMWF A subset of these are in Leeds on nfs portsdown_portsdown02 martyn ECMWF The files contain the spectral coefficients of the streamfunction velocity potential temperature and specific humidity on the original ECMWF model levels i e 60 levels as well as the surface pressure The table lists the periods that are currently available ERA 40 1 1 1957 to 31 12 2002 T42 Operational 1 11 1999 to 29 2 92 T42 8 2 UK Met Office UARS Analyses TOMCAT SLIMCAT can also be run using the gridpoint analyses of the U K Met Office produced for the UARS mission However this is not recommended for certain configurations of the model as the calculation of vertical motion from the divergence can lead to noisy fields These analyses are produced on 22 isobaric levels from 0 3hPa down to 1000hPa The analyses are produced on an Arakawa B grid so that the wind fields u and v are not given at the same location as temperature and geopotential TOMCAT SLIMCAT reads the data on this staggered grid and the fields are interpolated within the model This data is stored in hpcx home n02 n02 1rkd data UKMO and also nfs portsdown_portsdown02 martyn UKMO The table lists the periods that are currently available 21 10 1991 to 31 10 2003 2 5 lat x 3 75 lon 9 Running the TOMCAT SLIMCAT Model 9 1 Parameters The table lists some of the main pa
70. the points regularly spaced in longitude Instead of this note that we can use the Fourier decomposition of V V A u D Vm p em to write i Ai J Vdd XVm 11 f edd 40 Ai 1 A i 1 If we note Af 5 A i 1 Az the longitudes of the Gaussian grid and AA the longitude increment of this grid we can easily transform the integrals to A J EJA Lima oF 41 Ai 1 im 38 As E 2 A J T d sim EE eee Aii m 2 and finally we obtain for the meridional mass flux F a formula analogous to that giving V with a multiplication of the Fourier coefficients V by precalculated factors Em Van gt m 2 After calculating Vm as indicated above we therefore multiply them by a factor precalculated before per forming the FFT which leads exactly to the integrated value of the meridional mass flux at the interfaces F Y F e m This method requires the prior tabulation of the Legendre functions Pm 1 and their derivatives 1 12 dPam du at the latitudes of the specified box interfaces in the file LEGIyy see section xx 4 b calculation of the zonal flux We proceed differently for the calculation of the zonal mass flux If we consider the same grid as before we will first calculate the divergence of the zonal mass flux obtained from the total divergence of the mass flux minus the meridional divergence The total divergence of the mass flux is calculated by using a similar app
71. turated at interface N 2 so that RH gt 100 then the first criterion for cloud occurrence in the column has been met and this interface is set as cloud base If RH lt 100 the parcel is lifted to the next level tested for supersaturation and so on until cloud base is found At cloud base the lifted parcel now lies in a regime of moist adiabatic ascent Excess moisture is condensed and the specific humidity qp and temperature T are iteratively adjusted to their saturation RH 100 values qp Tp using dp sat Tp F dp p L dqeat Tp 1 ae L T T 9 z y from E AT Aq gt Aqz dq q Cp and dqsat l T AT q dsat dT 13 To achieve convergence towards saturation values in each subsequent iteration the new values of qy and T replace qp and Tp Finally cloud liquid water content is set equal to the total condensate The lifted parcel is now tested for buoyancy with respect to the environmental air If the virtual static energy of the parcel 54p is less than that of the environment Sve then the parcel is not buoyant and the model vertical column is set as non convecting Sup Sp Cp Tp 0 608 qp LWC Sve Se Cp Te 0 608 qe where Sp CpIp 8Z Se CpTe 82 If Sup gt Sve the parcel is buoyant and the second criterion for convective cloud occurrence within the column is met The final criterion is the horizontal convergence of moisture below cloud base This is calculated using
72. ucture of the model The cycle loop corresponds to the frequency of the forcing winds and temperatures typically 6 or 24 hours The iteration loop corresponds to the basic model timestep set by DTO e ADVEC Calls the subroutines to perform advection e CHIMIE Interface between the chemical model and the dynamical model e CALSUB Sets up convection and vertical diffusion terms 23 e CONVEC Applies convection and vertical diffusion e FINCYCL End of cycle Writes REST file e FINITER End of iteration Writes PDG DIAG ZON files e INICYCL Start of cycle period of forcing analyses e INIEXP Initialise experiment e INITER Start of iteration DTO timestep 9 5 Advection The choice of advection scheme is set using the parameter IADV PARAMETER TADV 3 SWITCH 23 TRAD Comment none MIDRAD Midrad CCMRAD CCM SLT ZH Morcrette 9 5 1 Prather Scheme LIMIT FALSE ADSLT 20 9 5 2 SLT Scheme A switch to ensure interpolated tracer volume mixing ratios are within the max min limits of the original values LIMIT FALSE ADSLT 20 The latitude at which the coordinate transformation occurs is set by DLATLIM 45 CONV BACKTRA 37 The vertical interpolation can be done in a cubic arg 1 or linear arg 0 by specifying xx in the call to POIDZ CALL POIDZ Z0 LGO ALFZM1 ALFZPO ALFZP1 ALFZP2 LPRO arg INTER 56 To preserve tracer sums DLATLIM 45 CONV xxx interpolation Recommended and default parameters for S
73. utines However for ease and to keep the UNICAT MPI version as similar as possible to the standard version the model will now only deal with rectangular patches of equal size In practice this means that the model grid LON LAT should be values which are convenient multiples of the factors of the number of processors requested Because standard grids tend to have 32 64 128 etc longitudes latitudes etc and CPUs often come in multiples of 16 in effect this is not a restriction at all 15 2 Defining the Grid The model resolution is set in the deck PARADI Additional parameters not used in the Scalar OpenMP model define the decomposition of the grid Parameter Meaning NLATMX Number of latitudes in a patch excluding halo NLONMX number of longitudes in a patch excluding halo NHALMX Max extent of halo in gridpoints NLATPP Number of processes per latitude NPROCMAX Maximum number of processes NPROCI Number of processes along longitude row NPROCK Number of processes along latitude row NLONMX and NLATMX are the maximum number of longitude and latitude points owned by any processor NPROCI defines the number of processors along a horizontal row and so NLONMX LON NPROCI NHALMX defines the maximum number of halo points that are required at a single latitude PARAMETER NLONMX 32 PARADI 68 NLATMX 4 PARADI 69 NHALMX 300 PARADI 70 NPROCI 4 PARADI 71 NPROCK 2 PARADI 72 NPROCMAX is the total number o
74. with the scheme within TOMCAT 36 Grid associated with a Gaussian grid A spectral truncation is associated with a grid in physical space a Gaussian grid where the physical and non linear terms are calculated When the CTM is coupled with winds from a spectral CGM it is natural to use this grid for the basis of the Prather scheme This will avoid using two grids and the problems of interpolation between the two The definition of such a grid is made by use of the following e The boxes are centred in longitude on the Gaussian grid i 5 I Ay 27 The grid is therefore regular in longitude like the Gaussian grid e The latitudes p are defined successively from dy 7 2 by using the Gaussian weights ww and the formula Hk Pk 1 Wk where up sin dx Remembering that Sw 2 29 k 1 K this defines a grid which satisfies the above conditions for which the surface area of the grid is 270 Sik T Wk 30 In TOMCAT the latitudes at the interfaces of the grid are given by DLAT2 Values of these are contained in the file grid data in home martyn SLIMCAT or the TOMCAT SLIMCAT web page for most common resolutions Calculation of the mass fluxes at the interfaces Suppose that the grid on the sphere has been defined with no further assumptions Now suppose that we know the horizontal wind field that will force the model in terms of spectral coefficients of streamfunction y and velocity potential y The rel
75. xing ratios as described by Prather 1986 This limiter can be switched on by setting the variable LIMIT to TRUE Note that the use of this limiter can destroy tracer correlations as the advection is no longer independent of the tracer distribution and so should be used with care LIMIT FALSE ADVEC 13 The following table lists some of the main model variables and the Fortran common decks in which they are stored 22 GRILLE defines interlevel pressure GRILLE defines interlevel pressure FOR3D centrelevel pressure FOR3D interlevel pressure FOR3D specific humidity at box centre MOMENTS total mass of box MOMENTS zero order moment MOMENTS fist order moment in x direction MOMENTS second order moment in x direction GRILLE surface area of grid cell FOR3D temperature at box centre FOR3D E W wind x direction at box centre FOR3D N S wind y direction at box centre FOR3D Pressure coordinate vertical velocity 9 3 Fortran Channels The models makes use of the following Fortran channels during execution Read file TRONyy Read file TRONxx Write PDG file Write DIAG file Write STAT file Write TRAJ file Write ZON file General I O channel Read file EVAPxx Writing reading convection diffusion matrix Read file LEGCyy Read file LEGIyy Read restart file Write restart file Trajectory restart file Read forcing files Read information about length of run etc 9 4 Subroutines The following diagram shows the str

Download Pdf Manuals

image

Related Search

Related Contents

NG 1, 2 and 4 For underground installation  

Copyright © All rights reserved.
Failed to retrieve file