Hi,
I was hoping to find out how often the UKCA_ADD_EMISS_MOD was called (eg. every chemical_timestep) by the model and how often the lower boundary conditions are reset (it is done within this module so I am guessing the answers to those will be the same!). I have so far been unable to find this in the documentation!
Thanks,
Hannah
Hi Hannah,
The ukca_add_emiss routine is called from within ukca_emiss_ctl which is called from ukca_main. The ukca_main subroutine is called every dynamical timestep (i.e. 20 minutes for UKESM-type climate configurations).
Within ukca_main some routines are called only when do_chemistry=.TRUE.. This includes the chemistry_ctl routines that contain the chemical solver used. However, the emissions routines are not within these sections and so are called every dynamical timestep.
Best wishes,
Luke
Hi Luke,
Thank you! I am wondering about the timeline of the surface concentration output as I have deposition switched on for an LBC species in order to calculate some lifetime characteristics post run.
As I can see the deposition in my surface concentration output, am I correct in thinking that the reset occurs, the deposition is applied and then the output is produced (and only after this does the reset occur again)?
Also, regarding producing AINITIAL files from a run, will this be done automatically or do I need to switch this on?
Thanks,
Hannah
Hi Hannah,
Deposition occurs on the chemical timesteps and will be applied during the call to ukca_chemistry_ctl. This is done after ukca_emiss_ctl so you should have something like
==================================
Dynamical timestep 1 (20 minutes)
----------------------------------
CALL ukca_emiss_ctl()
==================================
Dynamical timestep 2 (40 minutes)
----------------------------------
CALL ukca_emiss_ctl()
==================================
Dynamical timestep 3 (1 hour)
----------------------------------
CALL ukca_emiss_ctl()
CALL ukca_chemistry_ctl()
Do dry deposition
==================================
Dynamical timestep 4 (1 hr 20 min)
----------------------------------
CALL ukca_emiss_ctl()
==================================
As the model output is produced after the UKCA call if you output timestep-level information you’ll see a cycle in the fields where every our the concentrations drop as the deposition is applied.
I don’t know if we’ve ever run (or at least, not recently) with a lower boundary condition field also having deposition processes. To a certain extent, all these processes are rolled-up in the value that is set for the concentration of these species, and any effect of deposition will be reset the next timestep when the LBC is applied again.
Are you interested just in what the deposition would be? It should be possible to calculate this and output it as a diagnostic but then not actually apply it to the concentrations. Sorry if I’ve misunderstood here - it may need some careful thinking about.
The AINITIAL file should be an already existing file which may then be reconfigured to change/add fields etc. This is then output to ASTART which is then read by the atmosphere model as its initial conditions. Climate jobs will then produce “dump” files regularly, usually every month or every 3 months (depending on the configuration). These dump files are automatically picked-up by the model and used by the next cycle - you shouldn’t need to do anything here as the suite you’ll have copied probably has this all set up. You can use a dump file as an AINITIAL file for another run. These should also be archived regularly if you are using post-processing. They usually have a *a.da* suffix.
Best wishes,
Luke
Hi Luke,
Thank you that is really helpful. Yes I am just interested in what the deposition would be- rather than what I have now which is to not apply it to the concentrations! How would this be easiest done?
That makes sense regarding the AINITIAL files- I had previously created one with this job run for 2 years, however the most recent runs have all been quite short and must not have been long enough to produce the dump files.
Thanks,
Hannah
Hi Hannah,
It should be possible to calculate the deposition for a species and output this, but not actually apply the loss flux to the species. I haven’t done this, but it should be pretty straight-forward as the routines to produce the dry deposition don’t apply these to the species at that step, they are combined with the tendencies from the chemical solver and applied at the same time. You could calculate the deposition flux, save it to another array, then reset it back to zero so it isn’t applied to the LBC species, then output the array through STASH later.
I’m not sure I’d call this easy, but it shouldn’t be too difficult except for the STASH part. There are examples of how to output other diagnostics however - see ukca_chem_diags_mod.F90.
Best wishes,
Luke
Hi Luke,
Thank you, I have tried as you suggested and can get this to work with a 2D diagnostic for the deposition however, the other deposition velocities are given as 3D fields but with a 3D diagnostic, I am only able to get the surface 2D deposition repeated at every model_level.
I think this may be to do with the reshaping command I have used in ukca_chemistry_ctl.F90 here:
zdryrt2(:,: ) = 0.0e0
IF (L_ukca_intdd) THEN
! Interactive scheme extracts from levels in boundary layer
ddmask(: ) = (k <= nlev_with_ddep2(: ))
DO l=1,jpdd
IF (L_ukca_strattrop .AND. (speci(nldepd(l)) == 'H2 ')) THEN
WHERE (ddmask(: ))
h2dryrt(: ) = RESHAPE (zdryrt(:,:,l),(/theta_field_size/))
END WHERE
END IF
IF (L_ukca_strattrop .AND. (speci(nldepd(l)) == 'H2 ')) THEN
zdryrt(:,:,l) = 0.0
h2dryrt(:,:,k)= RESHAPE(h2dryrt(: ),(/row_length,rows/))
END IF
END DO
DO l=1, jpdd
WHERE (ddmask(: ))
zdryrt2(:,l) = &
RESHAPE(zdryrt(:,:,l),(/theta_field_size/))
END WHERE
END DO
There are also a few of these lines in the atmos_main and recon job errors:
exceptions: feenableexcept() mask [0x00000000] enabled. (mask [0x00000000] requested)
Is it sensible to try and get the 3D output working or is the 2D sufficient (even with the interactive dry deposition scheme on)?
Thanks,
Hannah
Hi Hannah,
The interactive dry deposition scheme only returns a 2D variable of deposition rates to be passed into UKCA. However, it puts these rates equally from the surface up to the level given in the array nlev_with_ddep, so the numbers should be repeated. As a 3D field though it shouldn’t be throughout the whole atmosphere though, just up to the top of the boundary layer.
The array nlev_with_ddep is output from the ukca_ddepctl routine, calculated in the ukca_ddcalc routine. How high up do the values go?
To a certain extent the only information that a 3D field will be giving you is the hight of the boundary layer as the deposition values are the same within it. If you wanted you could output a 2D field and the boundary layer height given by
REAL, ALLOCATABLE, SAVE, PUBLIC :: zbl(:,:)
! atmosphere_boundary_layer_thickness [CF] (m) {UM:00025}
Thanks,
Luke
Also, I’m a bit confused. Is the h2dryrt array the one you want, and is it used anywhere else?
What is the shape of ddmask? Is it over levels? You may want to have it as ddmask(l) and perhaps in an IF instead of a WHERE.
Thanks,
Luke
Hi Luke,
I think the value throughout the whole atmosphere might be coming from this line:
h2dryrt(:,:,k)= RESHAPE(h2dryrt(: ),(/row_length,rows/))
where I was trying to make the array 3D but it appears to just be repeating it over every k level so the values are being repeated to the top of the model in the output, but I hoped that including the WHERE ddmask would cause it to only go up to nlev_with_ddep2.
Yes the h2dryrt array is the one I want to use, it isn’t used anywhere else, just this line and then is placed into the stash and ukca_main1-ukca_main1.F90 in a similar way the other diagnostics in ukca_chem_diags_mod.F90 are.
Before I added the additional code for h2dryrt the module had this statement only:
IF (L_ukca_intdd) THEN
ddmask(: ) = (k <= nlev_with_ddep2(: ))
DO l=1, jpdd
WHERE (ddmask(: ))
zdryrt2(:,l) = &
RESHAPE(zdryrt(:,:,l),(/theta_field_size/))
END WHERE
END DO
Where nlev_with_ddep2 is in theta field size, I kept the where just to follow the example given originally. What impact would using ddmask(1) have?
Thanks,
Hannah
Ah - so ddmask is of size (theta_field_size), so this is basically saying is a point in the X-Y plane within the boundary layer or not. (It was ddmask(L) not (1) that I was asking about, but please ignore that).
The (k <= nlev_with_ddep2(: )) bit will return .TRUE. if the level number k is less than or equal to the maximum level of the boundary layer given by nlev_with_ddep2 (also now over (theta_field_size)). That should be equivalent for your h2dryrt array as it’s just duplicating things everywhere.
However, you also need to have an equivalent h2dryrt2 array of size (theta_field_size) and copy the logic completely. I can’t see an initialisation of h2dryrt2(:) = 0.0e0 and this is needed otherwise you’ll just be keeping values from the previous (lower) level where deposition is applied everywhere so the WHERE doesn’t actually make a difference as the values is isn’t changing to include the deposition values aren’t initialised to zero but are instead keeping the values they had previously though the loop which is why you’re getting values throughout the atmosphere.
Also, this line h2dryrt(:,:,k)= RESHAPE(h2dryrt(: ),(/row_length,rows/)) is quite problematic as you have defined h2dryrt as a 3D array on the left and a 1D array on the right.
I think you actually need at least 2 arrays, or potentially 3:
h2dryrt of shape (row_length,rows)h2dryrt2 of shape (theta_field_size)h2dryrt3 of shape (row_length,rows,model_levels)Which you then use in the following way:
h2dryrt(:,:) from zdryrt(:,:,l) when l is for H2.h2dryrt(:,:) to h2dryrt2(:) using WHERE (ddmask(: ))h2dryrt2(:) to h2dryrt3(:,:,k)
It is then h2dryrt3 which is then output through STASH. You could cut out h2dryrt and go straight to h2dryrt2 as you do already, but you need to ensure that this field is re-initialised to 0.0 every loop over k, and the final reshape is saved to a separate array (initialised to 0.0 at the start of the routine) for output purposes.
In chemistry_ctl the variables are only ever considered as X-Y slices which are then re-calculated every time you loop over k. You need to therefore construct the 3D field step by step. There may be a simpler way to do this, but this makes sense given the other arrays. It’s should also possible to go directly from zdryrt(:,:,l) (for H2) to h2dryrt3(:,:,k) from the values in nlev_with_ddep which is of shape (row_length,rows).
Thanks,
Luke
Hi Luke,
Thank you for your reply. With the modifications, there is a surface value for the deposition but every layer above that is zero. With the boundary layer output from stash code 025, I can see that the boundary layer reaches around 1700 m (15th model_level) in regions where deposition is non zero so my output is definitely not working as still only gives a surface value. Within atmos_main and recon job.err I have many of these exceptions:
[19] exceptions: feenableexcept() mask [0x00000000] enabled. (mask [0x00000000] requested)
Which I guess is what is setting the values above the surface to zero but I cant work out a way to get these to disappear.
This is how the code looks:
h2dryrt3(:)=0.0e0
zdryrt2(:,: ) = 0.0e0
IF (L_ukca_intdd) THEN
! Interactive scheme extracts from levels in boundary layer
ddmask(: ) = (k <= nlev_with_ddep2(:))
DO l=1,jpdd
IF (L_ukca_strattrop .AND. (speci(nldepd(l)) == 'H2 ')) THEN
h2dryrt2(:,: ) = zdryrt(:,:,l)
zdryrt(:,:,l) = 0.0
END IF
END DO
WHERE (ddmask(:))
h2dryrt3(: ) = RESHAPE(h2dryrt2(:,:),(/theta_field_size/))
END WHERE
h2dryrt(:,:,k) = RESHAPE(h2dryrt3(:),(/row_length,rows/))
I’m unsure where to go with this exception as I can’t find it being solved in any previous helpdesk tickets!
Thanks,
Hannah
Hi Luke,
Just to say this is working ok now although I still get a few of the exceptions.
I was just hoping to ask one general question about 2d and 3d diagnostics as the 2d appears to be converted from mol/m3/s to kg/m2/s but I was wondering how this calculation is done as my lowest level 3d output is three times the 2d output and I am unsure how this conversion has been completed!
Thanks,
Hannah
Hi Hannah,
Are the 2D ones emissions diagnostics? Emissions are in units of kg/m2/s, so I think that the output units of the diagnostics are as well.
The emissions are injected into the model every dynamical timestep (every 20minutes), but the chemistry is only run every chemical timestep (every hour).
Could you let me know specifically what the diagnostics are that you are outputting?
Best wishes,
Luke
Hi Luke,
The 2D diagnostic is just the surface level deposition and then the 3D is where is has been applied with the mask up to the boundary layer. I have used copydiag within the chem_diags_mod for the 2D output and copydiag_3d for the 3D output. This outputs the 3D diagnostic as mol/m3/s and the 2D as kg/m2/s. I can’t tell exactly where this conversion is occurring so wondered if this was being caused by copydiag being used as this seems to mostly be used only for surface emission diagnostics.
This is the code to create the diagnostics:
DO k=1,model_levels
h2dryrt3(:)=0.0e0
IF (L_ukca_intdd) THEN
ddmask(: ) = (k <= nlev_with_ddep2(:))
DO l=1,jpdd
IF (L_ukca_strattrop .AND. (speci(nldepd(l)) == 'H2 ')) THEN
h2dryrt2(:,: ) = zdryrt(:,:,l)
END IF
END DO
WHERE (ddmask(:))
h2dryrt3(: ) = RESHAPE(h2dryrt2(:,:),(/theta_field_size/))
END WHERE
h2dryrt(:,:,k) = RESHAPE(h2dryrt3(:),(/row_length,rows/))
END IF
END DO
So the 3D output is h2dryrt and the 2D is h2dryrt2 (I set the zdryrt for H2 later in the module). The surface level 3D and the 2D look almost identical bar the factor of 3 difference.
Thanks,
Hannah
Hi Hannah,
How have you determined that the output is in mol/m3/s and kg/m2/s? Is it from within Xconv, or have you worked this out by looking at the numbers in the fields?
If so Xconv will be taking this from the pp-field code and just interpreting this as the units, but this is arbitrary. If you changed the STASHmaster file entry it would give different units but the numbers remain the same.
Copydiag shouldn’t affect the numbers at all - it just outputs them. The way STASH processes them may affect things (i.e. meaning etc.).
Best wishes,
Luke
Hi Luke,
Yes I was reading from xconv - that makes sense if this isn’t correct!
I think in that case there must be a loss rate in s-1, following the line ‘loss rate [s-1]: zdryrt (lon, lat, species)’ in this guide to the dry dep : https://www.ukca.ac.uk/images/5/53/20180110_Folberth_UKCA_DryDep_V1.9.pdf as I am calling them both from zdryrt.
I think that means my code is wrong as it seems to be that these values are much too small and also the difference between the 2D and the 3D shouldn’t be there. If I multiply the output in s-1 with the concentration in mol/m3, I would expect the values to align with the standard deposition output but they are much too small.
When I look at the altitude profile for the 3D output, the diagnostic is not repeated until the boundary layer but instead decreases in magnitude with height.
Sorry for all of the questions, I am very puzzled as to why this loop isn’t producing the correct output!
This is the code for reference:
REAL, INTENT(INOUT) :: h2dryrt(row_length, rows,model_levels)
REAL, INTENT(INOUT) :: h2dryrt2(row_length,rows)
REAL :: h2dryrt3(theta_field_size)
IF (L_ukca_intdd) THEN
DO l=1,jpdd
IF (L_ukca_strattrop .AND. (speci(nldepd(l)) == 'H2 ')) THEN
h2dryrt2(:,: ) = zdryrt(:,:,l)
END IF
END DO
END IF
!Additional loop to output the 2D and 3D H2 lbc deposition velocity
DO k=1,model_levels
h2dryrt3(: )=0.0e0
IF (L_ukca_intdd) THEN
ddmask(: ) = (k <= nlev_with_ddep2(:))
WHERE (ddmask(:))
h2dryrt3(: ) = RESHAPE(h2dryrt2(:,:),(/theta_field_size/))
END WHERE
h2dryrt(:,:,k) = RESHAPE(h2dryrt3(:),(/row_length,rows/))
END IF
END DO
DO k=1,model_levels
DO l=1, jpdd
IF (L_ukca_strattrop .AND. (speci(nldepd(l)) == 'H2 ')) THEN
zdryrt(:,:,l) = 0.0
END IF
END DO
END DO
Thanks,
Hannah
I can’t spot anything that looks obviously wrong in your code above. If you instead output a species you can have an existing STASH deposition flux for (e.g. CO) and output it throughasad_flux_dat and through your method, do you get the same answers? Also look through asad_chem_flux_diags at the deposition routines and see if there are any unit conversions that need to be made.
Could you plot up the 3D and 2D diagnostics?
Thanks,
Luke
Hi Luke,
The outputs in s-1 are attached for 3D, 2D and then 3D at 25oN latitude to see the altitude profile. The 2D is less than 3D by a factor of three and the 3D altitude profile is not repeated with z but instead decreases with z. If I change the line WHERE ddmask(: ) to an IF (k <= 20) THEN statement, I get the 3D surface value repeated until level 20 and then 0 elsewhere (the 2D is still different) - k<=20 was just an arbitrary choice!
I tried this with ozone as this was already being output as a dry deposition (50021) and found that if I used my rate, multiplied by the output concentration (mass mixing ratioambient P ¶/molar mass (kg/mol)) the final rate in mol/m3/s was too small by a factor of around 1.510^6.
Having a quick flick through the asad_chem_flux_diags, there is a conversion factor mentioned:
REAL, PARAMETER :: convfac = 1.0e6/avogadro ! conversion factor to convert volume into cm^3 and result in mol
And the section with dry dep treated is as such when not using asad_columns: asad_chemdiags(i)%throughput(:,:,klevel)=
RESHAPE(dpd(:,asad_chemdiags(i)%location),
(/row_length,rows/))*
RESHAPE(y(:,asad_chemdiags(i)%location),
(/row_length,rows/))* volume(:,:,klevel)*convfac
Thanks,
Hannah
Hi Hannah,
Could you fcm commit your changes in your branch and suite and let me know the names of both? I’m going to need to dig through things a little bit I think.
Best wishes,
Luke
