From eee6f28173bb4179aecdf023262a54545e618cd0 Mon Sep 17 00:00:00 2001 From: Christine Johnson <74597224+cjohnson-pi@users.noreply.github.com> Date: Mon, 12 Jan 2026 14:34:48 +0000 Subject: [PATCH 1/8] changes from fcm for linear integration tests --- .../gungho/source/driver/gungho_step_mod.x90 | 11 +- .../linear/integration-test/nwp_gal9/ReadMe | 1 + .../integration-test/nwp_gal9/iodef.xml | 124 ++++ .../integration-test/nwp_gal9/nwp_gal9.f90 | 144 +++++ .../integration-test/nwp_gal9/nwp_gal9.py | 105 ++++ .../nwp_gal9/resources/mesh_C12_MG.nc | Bin 0 -> 128056 bytes .../resources/nwp_gal9_configuration.nml | 412 ++++++++++++ .../resources/runge_kutta_configuration.nml | 4 +- .../runge_kutta/runge_kutta.f90 | 49 +- .../resources/semi_implicit_configuration.nml | 2 +- .../semi_implicit/semi_implicit.f90 | 89 +-- .../semi_implicit/semi_implicit.py | 44 +- .../tl_test_advect_density_field_mod.x90 | 1 - .../tl_test_advect_theta_field_mod.x90 | 2 - .../tl_test_convergence_rate_check.f90 | 132 +++- .../tl_test/tl_test_driver_mod.f90 | 170 ++--- .../tl_test/tl_test_hydrostatic_mod.x90 | 2 - .../tl_test_kinetic_energy_gradient_mod.x90 | 1 - .../tl_test/tl_test_pressure_grad_bd_mod.x90 | 2 - .../tl_test_project_eos_pressure_mod.x90 | 1 - .../tl_test/tl_test_rhs_alg_mod.x90 | 109 ++-- .../tl_test/tl_test_rhs_project_eos_mod.x90 | 1 - .../tl_test/tl_test_rhs_sample_eos_mod.x90 | 1 - .../tl_test/tl_test_rk_alg_mod.x90 | 2 +- .../tl_test_sample_eos_pressure_mod.x90 | 1 - .../tl_test/tl_test_semi_imp_alg_mod.x90 | 88 +-- .../tl_test/tl_test_timesteps_alg_mod.x90 | 199 ++---- .../tl_test_timesteps_random_alg_mod.x90 | 595 ++++++++++++++++++ .../tl_test/tl_test_transport_control_mod.x90 | 168 +++-- .../tl_test/tl_test_vorticity_mod.x90 | 2 - .../plot_convergence/plot_convergence.py | 120 ++++ .../plot_convergence/plot_convergence.sh | 112 ++++ 32 files changed, 2095 insertions(+), 599 deletions(-) create mode 100644 science/linear/integration-test/nwp_gal9/ReadMe create mode 100644 science/linear/integration-test/nwp_gal9/iodef.xml create mode 100644 science/linear/integration-test/nwp_gal9/nwp_gal9.f90 create mode 100755 science/linear/integration-test/nwp_gal9/nwp_gal9.py create mode 100644 science/linear/integration-test/nwp_gal9/resources/mesh_C12_MG.nc create mode 100644 science/linear/integration-test/nwp_gal9/resources/nwp_gal9_configuration.nml create mode 100644 science/linear/integration-test/tl_test/tl_test_timesteps_random_alg_mod.x90 create mode 100644 science/linear/plot_convergence/plot_convergence.py create mode 100755 science/linear/plot_convergence/plot_convergence.sh diff --git a/science/gungho/source/driver/gungho_step_mod.x90 b/science/gungho/source/driver/gungho_step_mod.x90 index da004ab8..18732619 100644 --- a/science/gungho/source/driver/gungho_step_mod.x90 +++ b/science/gungho/source/driver/gungho_step_mod.x90 @@ -118,10 +118,13 @@ module gungho_step_mod call log_event( log_scratch_space, LOG_LEVEL_INFO ) temp_corr_io_value => get_io_value( modeldb%values, 'temperature_correction_io_value') - call modeldb%values%get_value( 'total_dry_mass', total_dry_mass ) - call modeldb%values%get_value( 'total_energy', total_energy ) - call modeldb%values%get_value( 'total_energy_previous', & - total_energy_previous ) + + if ( encorr_usage /= encorr_usage_none ) then + call modeldb%values%get_value( 'total_dry_mass', total_dry_mass ) + call modeldb%values%get_value( 'total_energy', total_energy ) + call modeldb%values%get_value( 'total_energy_previous', & + total_energy_previous ) + end if use_moisture = ( moisture_formulation /= moisture_formulation_dry ) diff --git a/science/linear/integration-test/nwp_gal9/ReadMe b/science/linear/integration-test/nwp_gal9/ReadMe new file mode 100644 index 00000000..911e9288 --- /dev/null +++ b/science/linear/integration-test/nwp_gal9/ReadMe @@ -0,0 +1 @@ +The nwp_gal9.nml matches the configuration from the default rose-app.conf, in linear_model. diff --git a/science/linear/integration-test/nwp_gal9/iodef.xml b/science/linear/integration-test/nwp_gal9/iodef.xml new file mode 100644 index 00000000..35ae0ee0 --- /dev/null +++ b/science/linear/integration-test/nwp_gal9/iodef.xml @@ -0,0 +1,124 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + performance + 1.0 + + + + true + 50 + true + + + + + diff --git a/science/linear/integration-test/nwp_gal9/nwp_gal9.f90 b/science/linear/integration-test/nwp_gal9/nwp_gal9.f90 new file mode 100644 index 00000000..2280e549 --- /dev/null +++ b/science/linear/integration-test/nwp_gal9/nwp_gal9.f90 @@ -0,0 +1,144 @@ +!----------------------------------------------------------------------------- +! (C) Crown copyright 2026 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!----------------------------------------------------------------------------- + +!>@brief The top level program for the tangent linear tests. +!>@details The default is to run all available tests - which +!! test whether the linear code is tangent linear to the +!! corresponding nonlinear code. +program nwp_gal9 + + use driver_collections_mod, only: init_collections, final_collections + use driver_time_mod, only: init_time, final_time + use driver_comm_mod, only: init_comm, final_comm + use driver_log_mod, only: init_logger, final_logger + use driver_config_mod, only: init_config, final_config + use driver_modeldb_mod, only: modeldb_type + use lfric_mpi_mod, only: global_mpi + use gungho_mod, only: gungho_required_namelists + use log_mod, only: log_event, & + LOG_LEVEL_ERROR, & + LOG_LEVEL_INFO + use linear_driver_mod, only: initialise, finalise + use tl_test_driver_mod, only: run_timesteps, & + run_transport_control, & + run_semi_imp_alg, & + run_rhs_alg + + implicit none + + ! Model run working data set + type(modeldb_type) :: modeldb + character(*), parameter :: application_name = 'nwp_gal9' + character(:), allocatable :: filename + + ! Variables used for parsing command line arguments + integer :: length, status, nargs + character(len=0) :: dummy + character(len=:), allocatable :: program_name, test_flag + + ! Flags which determine the tests that will be carried out + logical :: do_test_timesteps = .false. + logical :: do_test_transport_control = .false. + logical :: do_test_semi_imp_alg = .false. + logical :: do_test_rhs_alg = .false. + + ! Usage message to print + character(len=256) :: usage_message + + modeldb%mpi => global_mpi + + call modeldb%configuration%initialise( application_name, table_len=10 ) + call modeldb%values%initialise('values', 5) + + call log_event( 'TL testing running ...', LOG_LEVEL_INFO ) + + ! Create the depository, prognostics and diagnostics field collections + call modeldb%fields%add_empty_field_collection("depository", table_len = 100) + call modeldb%fields%add_empty_field_collection("prognostic_fields", & + table_len = 100) + call modeldb%fields%add_empty_field_collection("diagnostic_fields", & + table_len = 100) + call modeldb%fields%add_empty_field_collection("lbc_fields", & + table_len = 100) + call modeldb%fields%add_empty_field_collection("radiation_fields", & + table_len = 100) + call modeldb%fields%add_empty_field_collection("fd_fields", & + table_len = 100) + + + call modeldb%io_contexts%initialise(application_name, 100) + + ! Parse command line parameters + call get_command_argument( 0, dummy, length, status ) + allocate(character(length)::program_name) + call get_command_argument( 0, program_name, length, status ) + nargs = command_argument_count() + + ! Print out usage message if wrong number of arguments is specified + if (nargs /= 2) then + write(usage_message,*) "Usage: ",trim(program_name), & + " " // & + " test_XXX with XXX in { " // & + " timesteps, " // & + " transport_control, " // & + " semi_imp_alg, " // & + " rhs_alg, " // & + " } " + call log_event( trim(usage_message), LOG_LEVEL_ERROR ) + end if + + call get_command_argument( 1, dummy, length, status ) + allocate( character(length) :: filename ) + call get_command_argument( 1, filename, length, status ) + + call get_command_argument( 2, dummy, length, status ) + allocate(character(length)::test_flag) + call get_command_argument( 2, test_flag, length, status ) + + ! Choose test case depending on flag provided in the first command + ! line argument + select case (trim(test_flag)) + case ("test_timesteps") + do_test_timesteps = .true. + case ("test_transport_control") + do_test_transport_control = .true. + case ("test_semi_imp_alg") + do_test_semi_imp_alg = .true. + case ("test_rhs_alg") + do_test_rhs_alg = .true. + case default + call log_event( "Unknown test", LOG_LEVEL_ERROR ) + end select + + call init_comm( application_name, modeldb ) + call init_config( filename, gungho_required_namelists, & + modeldb%configuration ) + call init_logger( modeldb%mpi%get_comm(), application_name ) + call init_collections() + call init_time( modeldb ) + call initialise( application_name, modeldb ) + + if (do_test_timesteps) then + call run_timesteps(modeldb) + endif + if (do_test_transport_control) then + call run_transport_control(modeldb) + endif + if (do_test_rhs_alg) then + call run_rhs_alg(modeldb) + endif + if (do_test_semi_imp_alg) then + call run_semi_imp_alg(modeldb) + endif + + call finalise( application_name, modeldb ) + call final_time( modeldb ) + call final_collections() + call final_logger( application_name ) + call final_config() + call final_comm( modeldb ) + +end program nwp_gal9 diff --git a/science/linear/integration-test/nwp_gal9/nwp_gal9.py b/science/linear/integration-test/nwp_gal9/nwp_gal9.py new file mode 100755 index 00000000..aa141417 --- /dev/null +++ b/science/linear/integration-test/nwp_gal9/nwp_gal9.py @@ -0,0 +1,105 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +############################################################################## +# (c) Crown copyright 2026 Met Office. All rights reserved. +# The file LICENCE, distributed with this code, contains details of the terms +# under which the code may be used. +############################################################################## +''' +Run the linear model integration tests for the default (nwp_gal9) configuration + +''' +import os +import re +import sys + + +from testframework import LFRicLoggingTest, TestEngine, TestFailed + + +class TLTest(LFRicLoggingTest): + ''' + Run the linear model integration tests + ''' + + def __init__(self, flag): + self._flag = flag + if 'MPIEXEC_BROKEN' in os.environ: + TLTest.set_mpiexec_broken() + super(TLTest, self).__init__([sys.argv[1], + 'resources/nwp_gal9_configuration.nml', + 'test_' + self._flag], + processes=1, + name='tl_test.Log') + + def test(self, return_code, out, err): + ''' + Error messages if the test failed to run + ''' + if return_code != 0: + message = 'Test program failed with exit code: {code}' + raise TestFailed(message.format(code=return_code), + stdout=out, stderr=err, + log=self.getLFRicLoggingLog()) + + # "out" becomes self.getLFRicLoggingLog() when PE>1 + if not self.test_passed(out): + message = 'Test {} failed' + raise TestFailed(message.format(self._flag), + stdout=out, stderr=err, + log=self.getLFRicLoggingLog()) + + return 'TL test : '+self._flag + + def test_passed(self, out): + ''' + Examine the output to see if the validity test passed + ''' + success = False + pattern = re.compile(r'\s+test\s+.*?:\s*PASS\s*$') + for line in out.split("\n"): + match = pattern.search(line) + if match: + success = True + return success + +class tl_test_semi_imp_alg(TLTest): + ''' + Test the semi implicit timestep + ''' + def __init__(self): + flag = "semi_imp_alg" + super(tl_test_semi_imp_alg, self).__init__(flag) + + +class tl_test_rhs_alg(TLTest): + ''' + Test the right hand side forcing for the mixed solver + ''' + def __init__(self): + flag = "rhs_alg" + super(tl_test_rhs_alg, self).__init__(flag) + +class tl_test_transport_control(TLTest): + ''' + Test the transport + ''' + def __init__(self): + flag = "transport_control" + super(tl_test_transport_control, self).__init__(flag) + + +class tl_test_timesteps(TLTest): + ''' + Test running over multiple timesteps + ''' + def __init__(self): + flag = "timesteps" + super(tl_test_timesteps, self).__init__(flag) + +if __name__ == '__main__': + TestEngine.run(tl_test_rhs_alg()) + TestEngine.run(tl_test_transport_control()) + TestEngine.run(tl_test_semi_imp_alg()) + TestEngine.run(tl_test_timesteps()) + diff --git a/science/linear/integration-test/nwp_gal9/resources/mesh_C12_MG.nc b/science/linear/integration-test/nwp_gal9/resources/mesh_C12_MG.nc new file mode 100644 index 0000000000000000000000000000000000000000..ebe518907aea38df3c4cc5b3839febbef455ba53 GIT binary patch literal 128056 zcmdp+2YeOfx_4)0C!r(KtA-AW3aE$ziV7$yB1&%|K!9jSkOToM_TI6JV((oPJNDjt z$KHGI-&*gU#DI9<-s|~pIlr^~Kl6WwWIcQDgdJep0fUk#`@h;mE65BjEwYK7htDl@ z_MV(UWpiemP``BQ{L;di#Z$@)O3S8(+cT3_?2D#O4}IRU75iz0Q?mC@qHp-;$C+6; zt6+t*cBM1tlvET?pItn)prl8^l%kT7@`71KvkR&|^*la8)^qju&ab++<7z*ypDNrQ z9)pG}-e0+I^a{7EbXK`@eC=P)DtC)s;;iY4AG7-I;k@0SkISEXukdhWclUW&{k_@U zy|?3PKklEq`#g85?%A{2bFc~*UoUHR|DVrA*ERb1e>(X&R(Q;NR=xhB-hVx-d_E$d z&Xs;$)$`GdvlX6?o&_sBAD*wUSM}Fowa1}n!Jp5`n%)2B^RmL@@V6iT&&R>%rt*Jv z7}`P2%CC#`pf1@AYD}3{FuSZIe9`$Ds5i5yd`3Y<*{rgXvguicsOI#dvYAB{v%U15 zTII856wNN4QkZ@n)U0;9=Qa5{EH0Z`JjHvbl$Fh%T2MZ}oD2Tvd6Iqs)vmm++G7(B z54S~sp8feaRGu}QQdU}CF?-IGin7@S#ig_6RAfK6AA6+FYWAhruSd^5KIYVX{wR+3G)^H&**~;TQ8NTdikFap@l6*VQg9 zo9zoN{ZfkZ(#MsEqW&WngycVUO6Qc66cm@26qZi4Tj9RSo5)5*d(J83^<7+8QczM@ zQCu-6`x~kESH9M~z5dj)nT5rr1x53CT~_-#Xjb)j{`H0$rLzi4i%Nobn!}Q zi{Bsq?m6-Gk_+KAOlDx^?<4=}`CirQe6t~ZZR8j2RWy5kKATUSUCh=+vkOb|eJ=7V z%JPdV%JZl3CL@1J+3eDy+4;UE$~WXwM|H~ivaaAe!PKJgO}26;A17wNMf<-l4fxMW zJ5@dJYrfN}UsJ2vsr#Dm6g*sJ4R@Kiie0K+i!1Nq#nB=`Z!^Gl10r_Y$|b5m5vk>%=ksa-IrxTIucVac5AH@onh{qe_3U*THq zIa|2mp5d_%PdvAU*I)Q`%)VY%`8urnC4A#r>Faq_d-kmMocT+A<*{jg#oV&&D~_we z3SNx748w=Zz696wd3(ft)~fp24BywnYu4AQ{*`;w<27A6wUFOBSN#1TI9Gi)SoQB2 zt^9o?+>`!(;cpSuzFu%FKDHw-)UCi8@_3; zsc+3{-@I1!c(<%PU$ttqvTMCc*YN!|NLKf4l~48O;~QSzRgbUVv;K+4xAn@8ub(D= z3g3Tchu=h2cCK6LT@VNL!rQp;pEdg09umH`%lSUV z)8#MgQYq4vXRoBhul4ex3SWs{mgkq1=4an9X78?U+co^ov)b?1^%mgo-mbRt`#n$G z?f>HY@G0T@Zg?2`^;)C1;g!!t=Wf;Ck@sBnJMtbqyKNlal>b-WllwYd`8|1)N{#$# z&vs9rsa1ck{r~m;y!vC(qw+m^HN~_G-#h;9efo+YBm2Fv`kFn$AF&jKkMXbGzxlW| zTj6tuk5%>i^6K8RulQ@LeuJ*|{yY2PS^4+tzxll0T=8PyfdE@s;08=YPBEy{cXZmHSlY-|w~R z*V(Gxr~Dh<%UAAmxbL}E{F(?4dev*Qa*xXVU-*8$a@U1_eeS|@Q1!jT9sM&a-#`5A z@B3B1WX)ca|EBl!l{=m9^VevFuh$jd*ZXJZmAm{mzpu}G3Kt&TYkVG7;av6Jb4{FA z_ulg#dT(FpeA9h@aSmT6RqyY8jQ@$p`5$K^BR z=skYkGF)-KZ#e&5-S0nsoyy<5&wsQrXDwFvny7qZwc`8ydfc}b)Q8IF?w^03UnR0; zZ}j;KrApiX;QRbXeeBv)eLDYYs@g^Py8~|EfBv52?hRpW_@{Os_LrY=sy_c3?_TvQ z3QBse;r^LP?7*wxrRVr>^`CRnx6Wa||BFB8jL5!)@W;h|8?VvNIn^HBfBn}D(SPM< z8^11B{@JGb^U{O!RlU#o@A!Ep``2{*okO*UChgYv=b4__pIuh^7=M1R{g3=SlYXqK z_ay(3pJ#d&RJ}Q=_A^QLMY6h|XVTB>&DDNw*kQFl&!oS$>OJHCugnf{TqIs$xLDo_IdMk{=15Os$QEb@AJR#=b6kTc4XJLf8@U_?i!whs_z}{+Qd8I zZ(oyDzvRF5=b6kT95H==H2%*f8U@^QCwUIZtC&xOxROVHb0oc(u= zl|O?zw+6lKoZEojcFyZS<w?~P&iSCXopW2z+s?Ti=xyiR9`v?z?f`n*IbR0u zVAl~kf!^-o+!^$?bM69q+c|dyz3rUW1HJ8>*9X1roV!6j*g0fyJLkg-V^k;bDj-)+c}qm-geFvptqg#9MId&c`wl0&Ur5A zZRb1>^tN;U4cx(QKI{#8yNmNaptqg#zM!|A^M0VWo%8;nx1I9=ptqg#0?^ye`9RRy z&iNqF+s^r5(A&=W5YXF>Ij%*zgWaKU80hUT&WD5EcFsqD-geGMg5Gw{M}gjU&PRjZ zcFxDZ0r%zD!~XF4(K zP3N3x$*eb>bEZDC-gM5HgjsJo=j3~4z3H5jkD1+p&N+F5S#LV$70`XnB9TS zIk`2vL+BRKp93esiRpYUd2-r0FNTF+=X@TV0(Q>l!>M5Bd;y#WcFq^V>0sx45u5>b z&KJX(VCQ@ZoCS8ym%`a##~g0}ZPV^D@^WwocX7S~^tN-p5*CA<^Hret;&i^6)Z5Pa z8qnL$`C8E1&iOjf+s^rV(A&=W2GHBi`9^RDJLY)v><*z@LU$AB?JmwYgWh(|OF?ft z=Uc!X?3{T%gSVaYZJ@WiINuI>+d1C>dfPeQ33}T(-vxTxIWGfuuw#w~W_JkP-E{YW z-tOXjFX(OOd>`m-=X^itZRh*|=xyixAn0x9{1E7E=ln3}ZRh+5=xyixDClkH{1_|) zJLY&;c8AbCPWJ@p?Jmwwg5Gw{Pl4Wc&QF8hcFxa$-geHrDNcX56d^tN+;4fM8iejW6-bAAK#wsU?H^tN+e z4tm=;zXf{RIlm2h+d01jdfPd_3+`aY92aMI2;F;h?}Og%;`{;VZRh+U=xyix5$J8_ z{4wZl=lluiZRh+cEC)O1&p_|x>HImVx1IABptqg#m*5U|%<-J;4x#&s?rYH7U7WuG zz3rU81-!AbA3{8 zJLk1PZ#(A(ptqg#+Mu_cb3@SE&bbljZRgw=^tN+81Kh!`2{Z-0-Nm^X=xyiR9Q3wx zZUOFK=iCzXwsUR;dfPd-2EFZ^+koD7&g+2QcFya9-geIUptqg#Mc@u}$*pi3v)=6D zoZQZ=H=T2G2eaOE&dHt3deb>4cQNZt=bS8K)|<{bxtm#UI_KmbX1(d0lY5!egu=A`3bObxkdO|0#^K4`240fJv z0$sq)vrVBZ*m>3q)&o1wHiPxS&a>Xo4eUJgvAhoKJlg{F)O~a41A5wd))(}&^Q<4} zY3JFNpr@T@TY;W-p7jSk?L6BW^tAJA8_?6vvu#07JI@Ayo_3x+4DMj(*&xtU_kpk- z=xOKK_MoSoXM;gcJI{6iJ?%W(5%je4Y$wpu&a)w)r=4dzgPwMt4Fx^zJlh5IwDatF za0ffjb_G3k9|pUDo_3xM2R-dP8v%OSc{URCwDW8f=xOKKXwcKnvoWBjoo8b~Pdm@X zfu44rjR!sLJbMe=!OpXZpr`H=paAr=^K266Y3ErX=xOKKWYE*jvnimboo7=)Pdm?w zKu^K2i`)6TPfK~Fo+_5(fbJlh}iwDar$(9_Pd1)!%LXGt#U4tAa$1bXUz zARG*O+Ie;e=xOKKp`fRoXNQ5FcAgy$dfIt*1n6n!*^!{9oo7dZo_3xc4ens)*)gD} z9cRhf=?-?D9S3^qek>dhdfIt*0_bVy*@>X1oo6S3o_3y{40_smwh;8R^XwGR)6TO~ zK~Fo+P6Kza^Xzoc(~h&GO?HRSd3Gk~srwml7U*f`+1a3{oo9p z13m3LyBzei^Xv-H)6TOiK~Fo+t^z&nJi8k7wDas5a0ffjt_3~qI7>Fo?hrc8GNb7A z#NGE5nbFL8(s`B{!>lKrXPL3gdeV888ON+AooAWx%zDy!mYKk;C!J@RiOhP^d6p?) z)|1Y&%p_)apz|zK$m|Yuo@FLycd!d*cXD<;TnBf+U2p@~dA1C01Ut{}h9zL<**$O* z*m-s@+zfV}-3Lp-&a?aB7O?Z|0k{?HJbMss13S+ig4@B4v+#FZaoTzI2i#^u2zuIi_7do6=h@4kr=4f7fSz`qy$X8TdG;FUY3JGNpr@T@Z-Ab5 zp1ld~VCUI#(9@2ysAG1A(0TSY=&Ac#@DAu{=h?fUr=4f-fu44ry$^cYdG-P5Y3JF8 zpr@T@AAz2Bo_!2@+IjW~xPzT%pMsutoJBo=9YW{X=b)$VpTQTPr=4eCf}VDseFb{j zdGd6v{+b_d-# zOX`9<=+0SE5A@XCv-+T?oo8zSo}u%s0qAMRS@=7y!D;7NL(o(AwV@H{Y3EsE(9_Pd zCZMOCXH7v*JI|Vdo_3x!2R-dPYXN%NdDas2wDYVLxPzT%twB#a&Z5HX4x#gG9ne$v zHn1+}Y3Er!=xOI!ThP#fSz`qbp$=_JnIB{+IiL)+`-PXE}*9! zXVFYxhtPSp9_XoiS6CnPwDYVR=xOKK2B4>%XB&c^cAjkndfIu`9rU#GtOw|6=UGqC z)6TPv!5!>8+XVEq<1CsF><~K7dV!w0Zwi}%o_3z~20iUO+Z^JIDwy2+roI8On+#d#{|EzU)tw>VD&y~TMt=q=7OKyPs_2EE03chFm$_W-@cxdhxn znB%V59dwWK=w^c6;#>-Pi*p(1EzYw*Z*kre^cLsYptm@egWlp?0eXw`9MD^w_X54e zc`mqvFvnf8JLn$OpqmGJi}QTYTb%a>y~TMS&|94M1--?2KhRs8_XoYj`2f&coEL!J z;(Q?JEzSplI|y^!IlF`IQBAsoL2q$B1oRf?LqTtGJ`D60=fgp7aXtd{7Uv^DZ*e{f z^cLr%L2q$B2J{x^W5FGSIqsC*LHDQ@-Ep9|I3Evsi}MMfw>X~&dW-W(ptm@m40?<6 zLeN{BPXWEf`BczboKFM2#rbq_2VssoW_Qp%s!ewW=q=7?g5KhM7U(U`XM^72ya@Cb z=W{@BaXuIH7U#vFw>X~%dW-YB`6AF;oG%8w#rYD@ zTbwTiy~X)5&|91@2ffAl3ea1euLQlt`6|#`oUaCV5azgjb_d;~x^&lo-r{^M=q=9I zf!^YLJ?JgYH-O&ad?V;B&PzaValQ%k7U!EmZ*g7xciK5W0r#Yx^OJCI+BrW3_obcl({O*$h3Qgd=}inU7Vi- zz3rTz2fgi_UjV)BoL>aJ?VMi%z3rS|2EFZ^Uje=CoL>dK?VMi&z3rS|2Y0Yz&Wz6P z5V|+$-UPkf#d$gCZRh+J=xyixHt22V{0``C=lm|{ZRh+R=xyixKIm=d`~m1~=lmh) zZRh+ExPu*YW=wX6(0xqz3Fz%E&YyzbcFv!H-geHPgWh(|Ux40r&R>GwcFtdc-geGk zgWh(|-+O2#%6a2-S>1qfZp!n{3GaX=lm1sZRh+m=xyix3+Qd< z{43~f=lmP!ZRh+u=xyix2k33boQX)i>6|k$vpdi^XT||Lgk2J5ZVKM);+)ChT5oo7 z&g3%dP3N4+W7eC_Ia7mKZ#w5pO=i96oHMnU^`>*q)MnP3&N)+u*&XPdGj+io+{L*b z=cd)~x4Squ0KM&;*9N`qoEw7PcFv7JZ#(D4ptqfK6VThvxhd#v=iChR zwsUR{dfPd-0KM&;CxAQHwS-onx4SsE2EFZ^+koD7&g+2QcFya9-geIUptqfKThQCi zxgF?j=iDCjwsY>X zj-a=@IPU~{+c^&bz3rTL2EFZ^hl1XA&bxr#cFw~~>Uk!wS}G#iqQp+TBW$Y#(e&E{lF zXqsj#vJJFIldq#FA6lo`mTV8}rr7~iUTaZ1=m?#lGq{J@1-gQJn7pndUd!QU9h28_ z#A`RacQyT3yjCOkFnOIuyf!2EFna>8#mGHOUWd`9;2tKgyNK6X=p&Z=9tbjS-9_C&!7u>^~2lHW4ntPM` zfP0wx!hYZ$=KgR1xQDp_^m7mMAUGJ@!#o84Gf4Ibe;$1pc{sR-c?28@?qMDUM}vEq z$H1}R9_DdyJh+E>0-Ol$VV(phgL{|@;S_KW^Hewu+`~K_^jnbTndDjE9_HEjUq`Y> zv0t(;2!3sa2dFVc{yAG?qOaDSAlz&SHm^n z9_F>6pL>|s!wujb=8gD2OtMF`1a5+x!9C2Sa0@I>^H%aUa1Zl#xC7k7yc6yM_b`{i z-QXVPJ#a6$hj|~|5AI<;01twDm=A$|?qNOxkAi!ckKz9o$sW<;@B};w?qNO!PlJ1y z&%m?b9_DlKJh+GX0=x+BVZH<}!?HAAAzuaeFkge$!9C13;7xE3b2;ee9_HKd4!DQ; zF8)7}>=C^O@52Y+9_EMe5x9r>F?<5Qd;{)behc4$ zdzjzD58xi=kD%Z3G=C<40rxO}<+t6nN%n|-gWur~a1SyjV*kS)W}HBVYxXeX9LVLG zJp^m>f*e1>G$iXl$%^k>{VEZ(OkV9d|Gbj_b@A94!DQ87t96s zFz3O1a1V2D*azIh+!yq74|9Jw0Nle|fd45Zd-y&gJ_rs5_b?BEL%}`F!{Bgm5Az5( z65PW)3XTT%Fpq&_!9C35;COHk^8`2%+`~KxP6qcd7lMB7VV(-7fqR&z<9`*&9=?x= z&xEtUJ2(!}k&K^>72Yhj}9`0rxO(f}6oT%%yM(xQBTw+y?Gp-VS$w zdzg2^UEm((GPoPu!@LLX1@|!T1O42?d;lH<_b?yA|8;`%0*``wn2*8Z;2!1^ z@Fcj0`4l`2?qNOy&w_iH&%yKH9_9=1BDjb761)uVVZH*df_s>+fqw2`z5#E7dzj1d z|Bhr2-$%r6!#m&}=DYA7xQF>Zd;soYeh43ddzc@?C*U6Dr|=oLhxs{t0q$Xb315MG zm|w#;;2!3;pr3n~-@^~!9_Ekys929=58p?`Kf^EJ9_Fv`8@Px0JNyCeK_(Hz%dzfAD??hy4+7B)EtDQE)W4hy5{dEVzgLad14ihy4j~BDjbBNpLc_hy6nQFC*C_ z^rw=ifqR?+r^6ZG9`PV_XTu_J5BqcATyPKj#c&?DhyD3*0l0_#g>VtLhyBHH z3Al&-rJ$dC*k2A;fP2_qiT~XsdxZXK@)~fDtKeF=4&1~3dbk1H!~RBC0`6gd6Wk2$ zVZRh^0r#-K6>bCfu)iJd0Qa!J6Yc`{uwMrHxrhBda4)!r{eAeqOtMGlA0QtD_qZP( zf``F9>>q(g!9DCBgU7)=?4N)q!9DDsf~Ubf?4N;W!9DDsgXh6L>|cNv!9DC>0{z^> z{uOu?+{6Ag{J$dEBlK^OZ-RTg4$I*!a1Z;p;T>=f`*-0za1Z;R|pN`!7L1_ptvOz5(~J|CT>iu0^s(=)Wg_0QdL~euSUEJ?wvm zU%)-=e}&(`J?wvnKfpccGZDnh?m?ePfYLqeGdYmUHG9}+@}LIS>|vj&3Hq^zeWo_l z;hH_{Gj;KA%Qbu0XX?XRT(gJoDKZUUZLZnFKGP5yam^m~na0qBYxc0uG=*kdvxj}A zIkW)xpwF~~R^T4?t)UILhy6ORF1Uw%KHwMjux|(L!9DCd;J-P^9-;3hEdxL)NVZR0R z0r#-)i~lf^Jwm@FxfQraKj;rzgL~L-1KWao*bjh#;2!paU^{RR`|V*cxQG1?up_vK z{Z23h+{1on7z*xTzYFN+9`?J!Zr~pF!|^XB*(3BL$x+}QBVaU)0r#*U3**2&?8n0d za1Z;5Pyp^>KM4xKJ?tmL6mSpwsZa#&VLuI~gL~M|0R7y`L9Nfdc0_K2w*zW~%!9DEf!F+HJ`@LZwa1Z-^VLxyW`~Bena1Z+h zpr3o#9|Q-3d)Oa>{~08Eg#K{yNH`P@gQMUWI0EdCh2!C9us;D#g5$vcWH<#*1p8Cr zbXW-Xey+3NG_XG#&Ve(*{#-Z@7J>cwa3L%P`-|WbxB%=gjre9kUJUk^!xi8jm%){A z6}X4})o=~ChyAs19k_@6^>72Yhy9JP1l+^^Cb${g!+t5;0`6gdE8GU|VShW^0q$Xc zC)@?@VZRLjjwE}8{vPsPaF4s;KDZy;!~Owy5ZuH5A$S&(!~Pj~7Tm-BId~r2!~O+$5!}Q6CH(u5>=F7`$XCHVUWV7;b#M>+H{eZh z5Buft7PyD~+wcy!hyAX$p@IpHGA0S)P_1-vxo01a_T}ouGzyr zr#`I3HGA0SG=R0aW)J(EhR}#>_OQ=s3{ALZ5Br>^(2Q&Lu+M1@ExlU3mu>%xQBfw=nU>*-vzpYd)TiB>w|mPcY_VU zJ?uAxjlezZyF(9f5Br|5F}R2QCa@{EhkY;D4BW%MH~x!A_6YqJWFK&k&7m*!1NX4s z61D>OuZWbi~#qr9|@ztJ?uxr7;q2!u`mwY!+tzW0Qaz;2nFCC_LHCx z+{1n{Oab?>p9)3b9`@5&rk?axrgUExy zJr0CJ;81W6`@`UHa1Z+<;7D)}`=j7!a1Z-q;8<`E`{Uqva1Z+v;6!i_`;*{ga1Z;1 za0C2>FmrAsdvneHX3i3FbFO(z%$%FZEx5+R%(chL7I^)~M$wlcasB1z0G^bfu!E%z2qQLZ}UEKJ5q1+esX(KZ}S0iFsZlsAh`pH zcgTnEpUXAgQ4R88`W;EV%}2qAw)%Hs2;ElX{!)kW)y#&3DPEq~7LxWD%*i`93*~)Z6@koKEU(en`$B^)^2u zi%GrBkICIhz0FU^JxIJmeu{r7*LX(_$j|6YNWIO^$(f|y<`<;*>ur8XmXUg!Uy-v& zz0I#lA0NHVZ^+rC-sZPtIjOh#9qH>qZ}WR{4ym{K1GyK8cgP>{cOSf?M&wWQb4k6; zpUHWo-sUgld{S@oS8{JsZ}T^DA5w4gcXD4+Z}SgwKT>Zpmrt9!Ki7JbxiK8TwcccI z0t>j-o6OC?fn4Kl=H}q<^M-fSlxsdIx^;yMKyPzBxDY%u*N2NBKh19B#h|yj z0bBxln;XKVptrdZTn2iZ-QjZ3+w1{XfZk?LxDxa>H-@VK?~t3|Kb{%yh|d}O`L71O z$=Hv74d_kA{`zY{Z!+#hUI+dfW8Og^^cFAV$-6a0WuG!4Jhvai5!?T%t zFS#ApBeR)%AMb-oxgL$bkN5rb?&}`r1LR=RJLiEG%ye1z*E zq|wshHP3a}!+eSBv7~#LFO%a)_b^`} zeQmpk`6@YqbPw}2{AZKw5&QM>I{ie_JG6Lf%Ja&Fn=WHlI~&tM9w4K!~B_?Pr8Ts z3%NJx9_Fv)KBRk?zmfZr?qU8;?nkq<>V>p0o_aO5USirS= zka-z6kZboK^K#%IuHA#o%Y}ovb`LU-Pmp&A*X}{))qq2}b`LVICLG4Kdysjx;Bc

JBgL{~5;0$mN)5qjYa1V1`I1Ajv^f^Bp+{0`Oi@-h1c5n{3huI#^ z1@|yJ;J-h~96Y^?s z4|7wv2HeB!1=oUmn47_M;2vggxE|cY+#GHI_b|7B8^JxyKClGb!|aRyK_q*`JV*YS zO5RQ29;AQfl6Nz>2kD>5$_{**=gQZQfGDweOppzb_#tvQfJ-mdEEzR-R(QTB4DTZ zIkF?@th>ALx{q{rCwpG^k`wL@g3h|z^STety4!b$%Yf%0UQYG^oppB?UiZOScY9vw>th>8x0Xpk$-v@No-M%j@0iK8WQ?eiEth>8x z2|DX;zZK}LyM2GqS$F%bL1*3Vw*j4Xx8D|Y*4=&p=&ZZ_K+sut`$3?y?)KY(&br%g z4|fAQ#XpjRL1*3FWe3n%cl#ZoKj?116X>jae;5Kf>u$d@=&ZZ_P|#U-`&~e1-R*~g z&br(03OegYRyZvx@9M~zzAxD7Dy1UCr&{=o;QJ}N#_M<`P{$M`_bk^N|Ea8R2Ay@cp8`7TZa)=t z*4@4cbk^N|8tAOM{dCY-cl#Njv+nlAptJ7wyMxZU+wTFBKzI8Rcn{bqX-&=qoppDY zQqWm<`!dj3cl%kOv+nkLg3h|z&jy`!w=V~sb+@koopraL13K$&zZdANyZv0yS$F$+ z@D)tLzY94Zbk^Nn_6D7Gx8DbJ*4=(z&{=o;{Xl2k?e_e+c{r?38Rm9tt|^?kNXWbXTv7oc=_Q!$Fy4xQQI_qwK0^fIPk?fT8=eN}n%sT7NE}6lk&brfQ z`jI;8PM_&P>a06`rX#7d?(~@!q|Un2XX=qU>rS7ENu3t}eexZtv+nfCN2JcW(qoPFbRl^L*N0#(Z%FhTtAWRqf5xMxqdR+N0;K? zmg}duZVs1$&bqtH<)AY=h5ib1F|bqUuOxM5r_f(T>da1|znavUokD*NsWUr;{#sIJ zb_)Gy#?<`wMeLFdI_ ze+%gBPWHEg&hBJ?8|ds#_P2x1?qq)l=d&{=nP zxd(LC-Tq$CS$F&UKxf_U?+2ZAw|@Y1*4^IMozA-3`?}LvcY9xVI_qx#2|)ZoeFK*4_Rs&{=o;w?Sv!?cV{M!%p_^;y;0Or{VA(=&ZZDybn6-ZvO%3 zth@b(ptJ7wAA!!g+kXr?>u&!EEC=1~KLwq2Uk;yv&br%w4m#^@{{`S2cC!Bx|6*?)0&R&bqTpoRB)}P9JAT zopq;=b4Z?(}g3&{=o-cx}*Gclx*?=&U<^+z4=XC;PZD z9Lu#k9m3oMbk?0+;--LeaHo%(0nY9ec4u$d;=&ZYaKAZvWbOOH*Kh3PO?(C90L+Y$Meex`+v(E#4@*JtN&jWq( zJgKwp^vMgP&brelFOoXzPM^F)>a06`@-nHj?)1qkq|Un2C$Ew^>rS7%MqUK&v;nlE zdx~AoX5J9mbA2D#7B+$oTrVZLFYXQGUAf+u2Hm(GO>$q{8?NJeCzAX8eKaZUA=Eh&{LZrg#*x_aOCptG*tHyCu*)%$h;optrT9YJSZ zz0b!@XI;H-2*{?&L1%aJzFpv9a3}8@20G&!EhBdYopp7$-9Tquy>B?^ ztgH8p0G)O9zLB7_uHH8abk^1Td>!bltM`oooptrTv7ocA-Zu_(*46vQ!}H)y-Zues z#x=Yzn#ioPuI^RCy@Sy%6y4La-UedVCDuHIJxpMg7h-yF~x*YLh*FJ_%}b+@^ov##DZ z4|LYm`{skrx_aN*{?AL1$gP?-bBkSMNI&bk^1TP6M2StM{Fb|JtNGdEXhJv#vGaOwd_Z zcRLGo*46vY2Ay^FzD1z3uHJVJ=&Y;voeMha>V1nrXI;JTJkVKJ?>irK*46ti0GxxX z_g#p88`7P;?;_Ay*T!%$=&Y-|T>?7m>V21j&boTvWuUXJ-gi0ZtgH830XpmIeOH3c zx_aMLptG*tcQxp&tM^?4I0sknyB7bhq&sn6crT%|uK8SlNb0OByCvU}I_t`PnTXU` zSMJNyC3V)7`!dZ*opt5DOb1eDUAZsQoz!_Ta9^e`sk5%!m)V}wSy%4M3@338J8@rT z5>A_P?M~iz2V76L6WI;!gd4ctk-QG>f*ZNsf#kk;87$#?Fv)%K-Eb4veqQd2?}3}S z-j3wH_+D7b^&pb_;``tht_PCb7vB%Jay@|LzW4#Sjq7bm?u#FU+qvF`*{@vgU-5o-xHv-uHN?~=&Y;vJq0@J>U~dx z&boTvGoZ7s-uEm#1iE_Pb5Ru4BiSk3_dKbyu094YfX=!;3@?Jtx_aMBptG*t_cG|L ztM|PEI_v6vuY%6Hdf#iHv##FvI_Rvc_q_o+>*{@P0?xtJ`VOSmC8>&k9PhQv9za$l0e2k6JO zJF#2%y;x^mE4a=hb=H;Lk{YDWx^iDqlhj#P?n`QsI_t`PNo`VRcjCUJ4(RMo+?Uh^ zo!yE1l6s)CuH2W@2c320zGN-HIk<8|(tw)=bInfP=kLYA+2>&ab3@Qs*IZ}>I_t`A zNn_AiSME!ifX=#dU(ytC4zApnGy|N2EB7VML1$gvtp(_;tGl%XoptrTR)BMG^}g2p z!Vx4ph5Oo&I^!DN8?D2vv###8F6gYQ_vM4mx_Vz*&{^(>V54&XI;Io1L&-) z_jLrFb@jeZptG*t*BNjQuHM%L|3cE8yssvn0y1Lu?ptG*t*9~;m)%!L8 zoptrT4MAsJy>BDXSy%7t4m#`VeLX;DUA?a-=&Y;vZ45XESMS>d|CyvadEcg>Gp^yi zQ7>klb#=GRKxbXOuQ%wdtM_dVI_v6vTY%2GdS4&VSy%7t3p(rSef>aZUA=Eh&{Z~jGWyX>^>&ktZF{IAAa$jaNsk5%!ml;LstSk3rMw0!(owzUjUEQI~?!sv|sWGD>a`WBKt*#!o2y)@e=!^lBg-<<7}UCHgZzA4)$yOG;- zy(HTw!^y#1-=1gN7oFK5^geGovqR{8UUg=N&`%}P_`?qFXEI=h2?8Q>gtu%Cs$uW?+X9CA<4Syy+N4La*;Uk*C!YF`05>uNs- zbk@~=FVI<6`?;XAuJ-dlXI<^*gU-6z?+rTZYQGQQ9ConZ7yqs#u2C+zALy*ByX+4- z>uP@h=&Y;#0?=7k`vXB|UF{D7oprT87MrMl&br!P06Obx zeuP^H=&Y;#9iX$W_IHBLy4v3bI_qk` z40P7j{%+7&SNnTFXI<^@1)ReU_V?l6p2Rh(L*5TM>*_8KfX=$wKL|SOYX1=EtgHRQ zptG*_kATj)+CK_9>uUcP=&Y;#bwl-BVUI)>q;N_I@DQL`shqj zXI<%|vq+r74)oF4__rlbpnDyj;kpU=G`s=Na@{oBCvTF^aosH2C(Fs_xo)2Alefqh zxNec{leftixo(;5lXu9MxNep0lXuCNxo(~9llRD1xNei}llRG2xn3vRCm)cnalLM~ z4}ZpX5!d-q6dh0Z5$LSDyL=2fvs36lA(sO?h5l1gXLbtxXQa;T6#CCeo!Ke$UywSp zQ|P}Wb!Ml~e?{udPNDyr)R~<^{|)&euv6&2#r(IU$t0xCy3=R)1;JT&`b-X~v+nGY$t88xoj#LC>a06`rUt3A?(~_Oq|Un2 zXKIl;>rS7kP3o*WeWnhe8g`=3)W!c&uGuL%hUwwO>+ph~c>u#S9I_qxV7IfC#z8&bS zyM24WIk?+*!2fcRouZ@2j-a#d?$QZ#*4@4{=&ZYa7tmRE`>vp~?)K|}&br&L4?62^ z-wkxu-F^emS$F#lL1*3VHv*i)PWIjLzk+0^=t!~$=&ZZD^aP!Cx8E3a*4=&+&{=o; zO+jbf?R$aFy4!CCI_qxV8+6v)esj=Scl#|sXWi}l0M212`@Z;JNwQOP1lbRC*4ux^)bk^N|An2^S{UFd;cl+%C=dhFg z_V`~#vQu<8IT&=--CcG7oprb05p>qwekagbcl#lrv+njggU-6!4+Widx8DVH*4=&> z=&ZZ_uAsB-_PYVjVJG|H_+L%3Q*;s-ipZFZ>tUG<;b0SJGK``o!mP z66j8!@H~e@Px&UQ7@ zYv)GkY*#bAc5a-`cD2%L=O*cFS3A9SZko<^b<%6+X6bBKH@$Xlp3ZhBrq|Ak(zylc zW786xO&^n1>1_Jgv`#zI$D~c#nLZZlq@C$wux{F!?w$|MraQMyXVYEVrJd=H?bFV5 zw+?A%y3=K8*OAn_6F8gh&^et=cj%IKraN>^JJTK3OFPpY)=xWA?`~;l=EDZyY`VjS z>1?{gMrmidL-({Z-C;@E^&s``3C^ZFY@E)fJ8Y76raSN&51r`_d<}%obcfB-&UA;~ zX=m!adD@xoz}H~tOn2y$cBVV@O*_*a?oPXYq~2SCv*`|7rL*Y{{lVFEhpp4ubcb!y z&UA-u)6R5<0cmHt!@#sN-Cqrk&{yyQH1z4#U#UbcbEj&UA;}(#~{;;b~{O!^>$mg4BB?IGgS; zDxFPt7@c;e-ec0vbceBNXS&0#RoU;d=8>dKhmW zN)I=fhpLCf3Fb|ZvfmR0HA)XRT8+}fP3EEWFv&cW9wwWI(!q3ri9=Aq-nt%4e* zhbdO0^f1*tlpdy;hpLA}q2Tch#eT=nf@f)4J;cwpwY04s;%D`8+Ex##Cz~ften!b- zo;3NHD@EIQQhJF#>!R#;^w}MyhxmID?`T^+#NS(aRom(zWutl0N!^9&9rUmxmM%dVcXR6tj4?3wyEb^jW^4-sUNT!?=IV> zUSKue-L_4=&}zJUY@7N)tMO*rHuXbR%S)VQ+tgL9UL@t%JS?bDj?rSPQI5?L^UyhB zsh~zV7R#(gIR?wkL+O2mc_=+UVjfDbE6qdcag})}y*+9ks-6<{&3jDBem^d#QF?g7 zYLp(HG!JFJSDS~@gRYmPhtk7S=ArcPw0S5!JYyb856_y1(!+D+q4eMdHL4yG&8=Q5 zWxv-6YLp(*k^Ku-iOTJtX>@_lA`Hep67R^zfF|C_TJw9!d{;%tPtn9rI9nc-K6X9`*`qlpfx* z8l{K#%|q$o1M^UN_|QC5JtWRG?;|Pu{js1%>ERQrQF{2)Jd_?jGY_SQ&&@;W;S2Ln z_WMio(7j@xphoH8E2~j@_}V;_9=N)IQShtfka^H6#?#XOWAikpYh zLkaUxdN|cQlpac&htk7>;PLZ1_B;7=I@KsWB!5n)8l{)y&*@a7^pgBJoobX`l0Tq--!xnkcE~wjRMtoTlwk zqN<>-W;I@E+orB=HC`FprmkT%URm3wu4y%1IoqbLWi?)T+onF zeTLO|m28{3j@5XTZJWBT)p%8Go4TIWc#qgNHRp-f>I;rd1FKPv(V6C<9Giybp&XM& z=Aj&mv&=&|293=_>D`%!(sL8@PK&{ zR_3Af(Aqqd9@?0P(nDMGPJ#;eM@P^0uP&}x(((#=EZVUT$!Jq$JvrH3Krq4Y4+Jd_@WnTOKDaPv@lINLmw9?mfj zrH6gyohxO(M+j<^9?r8GrH7H`q4aRRc_=+xU>-^j7n+CC!$szy^l-6xC_Ria52c6E z=Ara3#ypfBelqV8Df@k?phoH8GOJN~7;7F%50{&V(!&+zq4aR2c_=+xWgbcoSDS~@ z!!_oi^l+_tC_Rib4^E@yIFvC2Q z9&R%arH9+iL+N3rc_=;HVIE2kcbbRN!z}Yq^^j8Byt}0A_uYaTrH6a0M(JU;c_=;H zYaU7u_nC*%!~N!=^f1Rflpf}qhtk75^H6%2Zyrhy515Clhm?xuEs(O`3k5Yw4-Z<6 z(!)dMq4coGJd_?DHV>tT#pa>(u*5u+9+sMi(!(xrI*w`=1J)#wXbOk|P z$)CNZnf*=iojvDjuF$}!k(9!l>on}^c#E9Rl}`l@*- zJ?=0MrMI2tq3VgB2`0T=Qug~bL5sF)mu-iP89^NnyrH41oL+Rly^H6$t+dPyW z_Lzs#!#n1o^zg2EC_U^o4^|n1|BCLGw_0_|`m>9={^pc_%lAe@aQc}&6(o4#5=1JLa zekPdovT9rVosvzcM(H6XyLnQ2Njcs;DZQkeV4jp-QgWClrI(bP=1G(Loy$B_)2HMX z)F{2=u^Oe9yyl_wlFvL;J@7NZq?cdHexE3)QFka;M*6gCf~mm=n& z^itG3l>I)*Jd_?zHV>tTV&QF=JlYLp&InupTEY38By zP|7@%9!i^s(nA^ZP7lOGC_U6O52c6t=Arb^z&w;5&NL6Dhlb{%^w7vWlpfA952c63=Arc9%tO^fmI=X2 z9W2#;r*;>rQF=)2W}cK@QoEWbrI*w$=1J)#wX=CrdP(hMo|IlvJDMk@m(&jCN$DlE zy?Ii4No{AIlwMNXnkP+B5ByvcFQtikNzwnChB~F0)oD`d=7PF~)h(@VC8cgHsFjzj z$yCkvh*GFEf0k4!`^CPnztr?WO@Gw%P0jJ3<~UJv{HQsu)ICDY_ok9{Pvsw{x|dL$ z(%b4jR`-=s_Y>6ptsY?YKq+;)pdJ)zz9*KfbxyNnRXs$gP8n+TFsp}4sm~VF=U9EN z)gz?T=LzbOq2_yW$@+ZdXH$KFP@Qt2)fZWPv6On0pdM}Y7^^RlQeP^lFAFu_15DOq zm7iVpezy+Wu?dBp0KRU2GH(-P-Y&?TDagD-ka?#dbCw|UEM)P>}hcAoC$X<|0An!-CAk zg3Kj?%%y_NWrED*g3J|y%tr*7D+QUW1euQtG9MFUJ}$_7LXi2SAak`)W)^x?CcpPX zCco!HCcoE1Ccno+Ccn2sCcmdcCcl?MCclS6Cck$>CckHsvPoo0NLHpnlovS3;fIQ%b!fKhqVfC9*>bC^-+g9%hHQ)bA zp?+7qBlcSTp4IP5sXq|ZA6os9)gMc#KM~ZQTK!q5`F>go^%vrE@uk)Kto}+${k5R} z#_Ii6ACOWX6x82Z{avW}{#*+658`|Aqt!oI{j-$%UxNAg4zPsP#f>ia5^dtX5}}QfC*`$6I}Z)j6cpIR$ktt8-i3PfDFvV{S~vN~N#-B_F@oYhUNZYrg2CaBY_Zf8snu5&J1(~%4nP&(xxu40bE6A)T$gD5OY#_)yQ;^wEkl9F(*;tV2 z1er|)nN0L0Fn@@+fD4 z4!IV|>p$VzI2up5_Q_@AH`i$|+gV%}AGF(?pgqO~?KLnt(DiX)yI;wm{q6k&-wm(- zWBGfTFOgyYv;2Mg|A+bi^zr$>X#f9@9G}1O`2XSl|6g?c|Bu{Xf8+7D{ZCJxhb(T; z%USrJ8~j{YH-@YT+e7CB9X2hge)zbc(RO$};cdJh?8EC(-m!d+$FbuRT`$M3pJVrz zZU`UNhlIEF!4KRU?vnwC?iYKW@%DV)U!}CzW+w+mP=OwRCg8E;=+w;`1 ze4jiZ$Bs{WI1b0IACAMZ`zu->!snwbzRiTrPhK-AtMBm2tJ#a}o8ddJ$a7bTcIWwS z`My~9=ovNQb*I1trI!0HN8=T|E7v6p3DnIz!G1qPQ*HC|Y;|=#r+}7RK z`)GKr+ROHNH|k&6SG{{6&hPm7BYWPQ_(PoAKgWalPkj3N!+ArRoK@iX*eZzKr@E z(qPGDm(G304{f_4Z~C!s z_IIEQ`+m-^ID% z_>Ol`{@`$YogH6i$Jg2Q;Ou&EcD*^f9-Lhdajsnt&aQ{U`QdtSc0IUo{r+tG_p1xn z?*SLC-!vDl-yh@LaQ$X?cKtfLe&hUb{Wi6{-(0xgQeC*;raQaeoZWA6uHA1g+;20T z-EYqBw>UrQJKWcQ|L4f?{o&~Q9a(Juv)ca8ILiIM{fPS?zVDoB$NiW1JleWp$Gy27 z_n+e9W9!DhANR0sI9oSdopIbPiGS>R3EwAAwd=gQUFZAkIzJHiW9vw~|F(|A=kZ9_ zdHDW&cDp~v+x^+V?$5z?e-?=6iSJK44z`}y{kbaR{h2&EIzUIO8k+py(XK1reqLbf z;Wf4%*0%L!A-nhxpCb@vgRxe{Spee#<@NZJmkp zZJqJabGRMnVS538DE#t=WhEx&XK-P6PD9@-8HK-Z`-CRO8>PJe3ehI}Wsb+y0w_=(l! z2e>n~op;Tm3%0vDU2kuc+V(hGceB`b*7)^keDwa1SH+dRZY%*#E&k%)h{g{aNS3{&^qve`CDvM%&Tr(fIKG zg>hau{=fKe{Kxxn{IdFRJg@TM_>T19_&M&m6aO4KY-f#MkH$yux9i7;>;GCGuGcSoxPI6Aa6LEk;rhPKhwHtB59{us z!`sRGj9K`_hN|caabG+Y}$}zm-1Rj~DxJ|4#7X{#_ZbyYY7N{+$?)kKS*2 z$JQTPpN_3R$JU<-RgbMd$JU?ZSJY$AKgXVbxDbv#|Ago3_<8Br^G{eGj(z?*_W4Vn z=f)+khGU<%sjeo!8-zYumor=Sv3_aP8LplCpW}LDznH z;@O`D_K$ymFxrk@kH#n8-!Y5rU*CBDqSxEM^8AS6gPOXI?|idj@x9l&PIIPIC@}M; z`1dKJ?dbJre0YB$%c~pr8_ge$@7!_2HD_GZ)QNvTGujTX>*qN(e|Z0Q zmgnN*7xfqIPxSt-d$Qa)VqaGm{k~_k9lain5APok=VkWyNA33~;{A_)-_(A8;@{5i z(d}=-{eFMz#>YS7{zli!A0OXCM>xK9-Jc!b==%G+uZQS*j;_DxdWp_&)L(S`qV4GQ zXnc5o!mi(b8SanhdXBEY!`Dl6et&&bp1xxS#&)_ecQTY_rn0YKfbg3rI5FE$FBc=@qHZKZ_zxV zpStmN7F}P_br5Zb``5+Sdw74p`2PK;*Po-y9~9pw(e)f%chUM6oyX{SN9$UYpY~9% z_4|t5;@u}{EzaKlv2W0J%lqYa4fVA)Klsuc<1hB*CQqH5a?LG1PoA15+zP&iqX`d-vh@muflF<9`Uu-<+`VFhl#%2^;2?w(&8hZCw1)A>+kx(1^vba zf3FMUCb;mvWiI60@1oqaIS)VcTkmDw?Y;M<4u!7v4VvW6-l)erzLwv5#=1kReYxwh z>^uIE(|w*q|BbsJn&5(7l+W>5W5=iVUq3#-JHK=6{I1Awem~4`eh=9B zegDtT?;l>DT{2vs-7;LC_hq;~PqOQ?Q-aO!?Og9pUoU$9=?`4rniFRIx__hcb7h(G z^b$AVz?Ely*lv+aKX+Ku=R2)*t@_=6^3XDST)R{D@5)~Pa@TRs#kajXs)y^mJx87` z7q3?S*>^s@_BG|dG_2o6>y>}{-Sek5RQ{}Unw?eA^{w&U{yI}yxc>QD*2&jmha0f} z#zI#OY3kC?*Y<&Q77eaBqaKi`sHei|{-4fy7zt4~brsQe4> zY5!F-*J@z*rdbwtaqUhWeEK7w-{3mFd-+v^x=mC5ODm_k&dPuH{Q;l7sr>KqZ)xzf z@^gGuXG1gPmn=TFa9QP_-u|NJmMOnklkMG~P=1Gx%jbAT`IVMm``UTR-}CvcT_05b z-Y)yk*XP!EKit!@W4iLco1JCDS<3(Eii&$`DgW1vZ?(8U`6>0wte&p?Y!hB9wN?2! zhWyrGt@3kE@7Jr7^7FoSW9$CPKYl{L#f3BEH~TO1JN!ZZ|2#h9|Ni5XsqUqB!@+oc7}WiSy?gqE0 z_SNuQ+1!vi6W2dgwyPUj?u%?gw@-J&iuBAraCZ?mJo~0I-?$@(k2Ls2eHb^x+xxub z_$YT!kGr2ab9o2luN$2Ejm66U>9;56KJ12;Z*=-~*^4XxrpkHdPEr2rx20|9qkOct z5B^kdIZ8!T&<}jN9YweLmz2@KNrdp3{b|%>Rt??^=85pgblvziTa^Fo z$XSJ}DgX6To7_A@`3FiIzkinUgPN~=##L87?|VV{$sFZ}+(A8ebsAT3rt){+^5qqC zl>g4aqTS0V|Kr{*=TuPs*ZZrF`BM2m-}K?V1C$>$`1_R4xQfc>ePxwT&Rpe(+(A8m zxMF^{&y@dT!KeFnQvOeMJMVa1`9E*VGWdSw|ME=1ukKR*p%2cVUsw5P9dEw!88=7y zyl=Pi$?2^8kQ?cLmw*3%$p1TzPlf+*eE#S8?WyCFJWo14|MUD#*YWwQ=l37JKL3&X zYsumJ>u*^frrY|^&DMwP(fZ)+ImL(PmV@4&W4t}rcze$A;kjp@x91=qo{Oq`drtE9 z+~mV^R52f(t44c!PVwQnWwG*;&oSPfYrH+@D4*w^_R3E_2YGug^7fqM!*kPOAD*L@ z`S4s--`jJF56>-WK0L?l^7dS#e4cYEDWB(_M#@h<2Pr@KT;%OJ$%p5rUf!OgyggSb zKlz;E!*fe#Z_hExPd?WupXZ#nmCtj}cIERNR8#pp7nN0h@;OQ0XX3f3k`K>OgS|ah zDL?s~qI{lPo>M;0F{_oIe6CS`@;OKOJolt2KlvP_{N!_y@{`X=K0G(g@ZmXXg16@? z$>FC!dp)pL}joKF?9jl+Sb3 z!ha|~?SJGS&G9LBgyZwpQ5>KD=KPk)aDFTN!TC)-M`yS`cO2#Q`MUCXj$d+|h@ z{{9M|WB%g$P|end(izr=lWlz%W8W9ozaQ4)(|r@0f9yMS^v&iq{=|1a{o8WIAJ^}r zx7)a>OZzW;_w&E)*RyjO-*ZOJ`4267+xK2ozuwZ71%2O_a=!6?>!*DG-KV}&ZEr3= zaOcgF8t%y*{l2wtx~RyBm#;ZF{{DNLOP)&F$(;eTIj*P1g2O~1l7 zO22mOQ{zteP3tY~H~Qt5<9^z%?lq?JJ(v5=iL|#mj<4&x-BajCf7fK+Yv||M8xOz4 z_dWDfnH5)U_5(|I-SBy>w0M69?ce#(2RUE#gLizB@6)M2`7>%%KYLQSuY7}pyO*!u z>wS|J58k|U=#BRK`@i}2+jn2xxXupW_4&4K-}|Vt?|tT@kCr%jryo$Xz}oboGvnOB z)p``q@$ox;$Z0oxo_+kMe#jd?7s&HgDPQrfEgR0Cw$Rs@HT~{&=M42{{QTvNwFk@i zv(h_ey>3RIxS!77Jg{Vc;d^}FynVao`D{ph90srY{?v|JAM-=}yz?JidygMFa^Fi2 zt~}_6emrEzqH^E+{3X{kn7r>%#H*&V}>)rL*(v?EG5pZ_duI3+K0z zv-9ih{2rbsoL^_>*V*}Xc6~a#KI7bQeU5c@eLA~7on4>Ku1{yzXFQKxpU$pN7w)f< z;@oh5RdaTKIlI4{-Cr);UynMwzntA)&hD?n^Ca)D%Ma%!*M|?Ctq;!D2WRVpv-QD+ z_2DgN>w~lP!P)xYY<-C53E#KYzi*bz|LAYe569;e7ta&cpFcT1tB!npqJQtqj!!th zwGW@)@Ea=5)}Nz1zxsYplpEF`JHJ2L^{Ic)BdkAva(!0L^z~`$kKJEG?fxoj_gBC8 zxY_z+_t(&PfB*RYIvC$yw*J`q(Aw6Afwn$8Z0p1A@jSNvxZkf2Gi`l%KI8ge>)RSz z-_mS->uKxT0$bm1^S1u@=-&gg{0!^chK%c*t|U4F`}<2*9-QjZW_3EP$R%S;FEM|5JTAQN@i-^SZSYvm zgH;mWx<&;@b=Xxh-JNyd(~l1>-{V}R-rcKQS-`pbz8m&knHJ9N-hJQOGp;th#{8D? zxafURPLx}F;LPr`SJrcNI*(j-!}Z^|y6gvBd@4l(K^0)@4 zuStJ*>7-bLpEn*Cy)WeKcOkdFtFr7y|G+sLT=l69t|+w71=WE^MuqOTPHmYJ>PieyTpmET%qGm@Q<%^<;GMyrQpuX4r|h1Z{zaE?~8Jx+-$Qu<-GRV zv)u8&y#C^>S)*N!uNOS9DdlFDYvkHFPd>WE<-VZuHLdc0>vDhnUa3T-W{1`MXk2*T zEtWINg_V*_u+l(e8}JDLw>`p^qFs^tIK8<1^2P z5NYPtAFV^SE=B89=GLFL;&tWkU4OpL(2%8y7t%BAN^eT57(dhwvPR8tv~DHb>i<`f4+%x z|994(wg2V%^WIUeKO-|d|Ga1G%(TPLJ<)SeSl{gV=RJGhk)D4JW_bR2GsF7R>hN<< z@;NAa{?T)hJtzI~^UpySJ*R~88?8SD<8{TpzSwimH@O;uZ zUMIry&&dC!=byRuT$9(2J?5yY0Cq&WhHZNTbhX|8f2K_6Y0Gkv@O@)%EB3 z_;b_mKYwL*{b^@8TjS4l;p_W)Q?SbIt$G z=dZPg*PWwUf42U~=dZuI{reFg zYsTTvUCHMYeGW@LpG^PL^{3f6U8kI&3op`A2z$7TjKem zzM|t7a`bz38Lo%l=bxJ)e@dL2>3qBXqU+?3kI!fE{zm^`Vs_r+>^ySbz$c%9N)3={q+0wA>G^h?c9 zxzYYcxpscNT}R%oS8wMv*2uT>>!bIDoPKs3-m(5JjrSwUwd>Q{`s3}oJgnjR^wGFz zo~Z9|d@qfUZ*&}@{fcsJ{qc5RnNIPMZ}*4wInj>imN+*$Z_)9MjzhG+QLe2&-qweB zzGz(ZKD&MoUk70wiqBg#E_z?cpAzS0I^V9dKR!NH|c7Cd6I6j->;}g!q zo>a2!g1Sd=cj3`kspl<@9S?ltL^@(`fvAfbe~7Z?daBr=s3KSVSPxC|2vA& ze@}FH`1y78-*0rb?V|r)XNGP4ur6VwNFmiYP73@@G0@+iOOn~~+4Vj+-Z()D`g?Ur za!MhWk=gP3`*v9%kC9gj`HcKhIMK-Lc>P_xEKtzU-^I%U`nz~Zib$cTk=gP3J9=5* zWJ7;XFAM1J=_M&H1^!EiklFG2`+HfSq;Z-QN*Sf4P{zpY__BH*lrzdpp@LCS3i`W! zNisXWvfc+(jH*(oW>l9#4I{JTYwCSa%Q#&MwGI7!z$~D@511sgOIiO=q!aUMpr3xGcr5AyWR&qjGj{HW%QOpA0xBl z`|5qr&*(3O0meWnq#KzXKS=L`!Nw3N3^j&HVYrdm@n`FOaE@`V6h;{5Nnxas+41M= zeQ<$sp%gAME|$V5BeUa2>wPfBxI_w<8kb37tdZIAm+O6Sg>j`6t}?Ec!Zk)_$6u@W z!8qeODO_)im%wR#GajO)j7*nM%&B*Nd>3Sc` zFm98=?Z!+g++k#P{GECq%rfqh!rjI_QkZRIcKp41AKYi$FNHbATq(>mGCO|0-UkmD z3#72ncu)!t8JQiwNbiG(jm1(}Vl0)yG9$C&m+O77!gxdqD~(lBc+|-3_{a1nH|4Y?}K&5dMP|_Y>>i6BeUZ->3y)-*dm1& zj2ETwl9Ac*TlGHJW^9+j%f>5Gc-6@4_#Jv5>@;>s;WguRDeN{fHGW+9{FPnrJzk`U zgw-cVQ_Z8!A&xVTI;Y5L9(68}%{=PLqO5tjrFjH%Fc)=RK`tJ3K0z)Xb$&rE9`%WW zTs-QUf;sRCh=PJ#=Aten$i<^BEXc*9E+WXqqb@4Q#iOn(m;>)5ak3znxu}Z?a`C87 z5#-`g7Z>E>QI`s21i5(Br3JZo)MW&@c+|}VbKq?g zT>m8nb5V1@kz3O0r={fLQ9mxo#iL#!$i<_6STF}3=b?%yCphQSRYiG$M_o-+5O~zp zMMZ%}T|-n7c+?$5Yk|jks3pi{F6z?-xp>sI1-W?CX9#lfsOt!F@u<5C=D_1T)Dz?~ z7j=C>E*^CQK`tKknSxwA>V|?`JnFuJIq*0SX9;qdi@LEO7mwNra`C8}2y*eLn+kIA zs0Rt=z~elm338c>y15`1kGh3$0*|_-AlF&lN=hys_1S_s@Hp3t1-Z;cy+n|UN4-># zi$}doFb5v>azQR0^-jSYc%4KW!TIK#bQWy|9Lw$7{U=Ie5I*TabgtYkdScc)WJGU=BQ9>nF$|n{(1%kb}o-0|YsE zyf#pfgU4&>f*d?vyH+p- zwcUaoJYIW4kb}o-ZwltXL@wyj&tK6nMOLsklhs@!DnLVu8nNW5p4I|()HB4z2K8;?B7^#NaiKvyQ(RzBcM#_b zvQxz!f?S}!Q;-YPvjn+7eU~5?sP7i!0(E=A9AqCS?h)hy^=v^dP~R)a1?u|*xj=ot zAQ!0H3FaU>tC%Cm1?stiT%evO$OY>8f?S|}K#&X6Z3T0XolPtda}4T*f?S|}P>>7M z4+(OCdXXR(sM`qUAiJz6C&&fr@`7BTt{}(->WYF~pspmy1?tvr zFbCe#;u%3Mb5TDl$i<_6PLPX7?FG4b)N2K~c+{^8=D=Gg)(di(i~4y%E*|v;K`tKk zM&Sh>^(H~Cw|cjfIq)`%ErML;qJBY;i%0#UAQzANB|$D8^;W?gc+_tQ=D^z~SSQG3 zF6x&Bxp>sC2y*eLUlruyQST7s;!(dTm;*0Gq>8K}VVX^Pyg1JE1Zhr@-87dpkH}$~ zSDIhsHa$_=N#qj+L_twVFb`T-6cNmW78NH6=0Q&u#RT)9r->xTC@q)=EhEYb=0VGe@`8EL3ZkN59<-7;UocN)QAJc0B~7bIs|)5qYlxbHdC*$o zbiq7mZE=QR9<+|QTrf{vQBTwt%!4)%XA0&)8;VAPdC;>&W5GPA6HNs3piRXD!92}G znrJSV2W=r*3g$suiPnO7&^Ds2U>>xcXfK>;2kCUdJe=pwqO0g=+D+O+bTREI?Jc^S z_L24zy-fQ{2a3L?>C*ee05M1m7DEK{phLwl!93`2akgL{^c-=nU>4Niofw)jG4|jd+l*NgFjdC(igCc!)t#6)qUU>@`)F-b5FI$7K-m zPYCA0e^RU#%!9v1{3@6S|7r1zV4kPMv*J0yJosL$70iRbPOKNqga5qPAeaY#qn6IPc$x1L7-ze^7iU_6z*)#gF1!f&Y{Em-s>8 z{~}Vv&jSCD_)Rd6&U-?V;2+E*KcTUSRMnaXKXIJMs%_1KpU5V%Yg_Z+CyrNLKxiKP zL=KTt+nR^-p2#I~3t9MyJR+}P9{hYFzhEBx6NU1U^WYcM8&8#Lo}^z`T0}5UAyHJE zB$x;PWKm2o5B@2lxL_Xq5`sMD!7r)&%2Lgf^h-%g3+6dZlo4eG^Wc{gK`zQZrwipKHBGsTR4bR~vP5$+NRW%R5Q7D|XiG6f zG&OA{9V*C0TZ>^rxk+uK+=Z%@%X40$tr#xIMcaw91-WQ@agHDt?I6w-?O^w5vE@kc)N`7YK6E?&3nB+@$tU?i|(1O+JS*?jk|1 zlxLH~#e!TZ&nAgcf?O%jCW+AkU&^ydVvJBOntV37SIYieB5qSXBWU7M>FwGUXyP(4 zQ`>h0O^lVQ=IbPyxLn++?Yn~}u28*K%oZGvD+Obj2fa$%C71`jTJT=xL9Y?_2oNd|W5gOL89cdU2m%9(268Uoa1PgP0?j2c0113g$s4D)&R_I4Re| zjbff)9`q)`xn~}9l6XKc4?0;a5I30KEL|v=2fanPA4{31h4fbOpkN+!ig-vc4?0yW z63l~66AugKL8ptwf_YG#N51YP=Sh>+6H5g1pmoGj!8~Yfu}m-zT1zY!%!AetD+KeP z)s+9alzB#p>>@?z_2j~TyfjsCP2is(juSdJNk4})tI)Aa`Z=Z9gnCQ*xuj)<=1x8% z=9ZF~5X_ZFkV!UvUO^^v;^z}&GADk1K_+wJpC~E|=F~Ni{CY+v+00c?kV!UvAweeD z_=N?TWaAeRWRi_vRMZsAscRtlnUGAf3y6~inPlS^6J(N&e~KWJZ2aPaOtSGyh`NF~ zCyP@BnPi_NN(wT`#y?GvNj82dK_=Pwr3IN}>Z2Z#&nPlVF7G#o*e}*8FZ2UT+qhL*GRekoEXX7q-w86w#&06XBpbh} z=qs3$bHKSElWgWn6J(N&-&{CBHhv31CfQE36l9W(-%1P;%*lCwOpr-7b3HD|Bpd$; zK_=PwPYN=fz+WxMBpZK?I9o6$=b?>wRB%2ziMC>y;2P*G+KGn*cv^FSt9oR7hROtP4Jh#-?J#tju@lEt`T zB3+QhxZ&bj!JLfaJdjBi=i?kfCRxmVt{{^v#*GkUlEt|51es(pZlstX(gou<*JP5# zHQ)uAWHI+fK_*#@dr6Q<7UNzOWRk_W-Qs4!oQ&gKTp;ccoRcx)LNQZt9b6(V61NM+ zT`DdXw+Y5wCPs-Ff^lO-LQEHo<6MwQ7U$#&K_*$ueWf6iEXG|W$Rvw#R|_)9V%#+% zyI@YnaW2Rti*qtgkVzJEUnj^Ui*eTrGRb1xctIvvjJrYP7R<>w&c!%ER=Sud$Rvxo zZxm#bHBQ_l$Rvw#lLVP$F>bObAefVJoC`9^;+)(f$Q&n_`&L0FS&W+^$Rvw#Qw5o1 zF>abTNiZklINxNF#dYwyAd@WSeoc@`7UOmaGRb1xPC+JFjN2hj70k&v&O<`nDYynW zPn`Q%;x@r~#J@}2E;!HlcZ->V>j3{Aafdix;LjEv1X(HKUO^^V%ypk2lPvuE1({^w z&k3NKC* zTLqb9GuJjjCfWGg1({^yzbwck8~+tSCfWF}iV=c2ONt$WOtP74ry!GT{9S@fvhiON zWRi{lx*(Hm{N3U_!JMaxX@X3$nQN*blWhDcf=sgUZxv*cjemJzNSJI=PLb6Aa+)wZc~T8)>@wyASjjaSyTsViHZTgtxV5!5LAkk@LI zeaUAY%0A>b4`uEX%|n^9rg;UVD07rH4`q%r=Aq2d z%)D(Z(@bRj_U9 zYF6V_v~B9@R^wH&ZR#3UW)@(U9c~;1U1S&oo+SCzSK4kWgpHk4`uE;=Aq2l z-8`-fa_b3dlsW2KjWS0A^HAnE(>#lsQ^ijWWmC=5bw+yI4@8%(293lsT4~ zhmyO@Jd`<>n};&TPV+eDZ3O3=x|7v-ZEc&nv(-GQRWz8HOd@A%|n@En0Y9u zu6ZbPj4%&nj`PeznPa4RD7iPA$92IRy971L9Isi8GRN!Yq0F({Jd`=!Fb`#pH_hXm zUm!T&)MKp1yU@0&FR>c$BHO0E)M~toZJYWstMNwJHuYGm@e;O8o!x4#3-;v-L5;Fc zS6Yp-FISm|vJY3Ahcfpy=Aq1)+dQrda>og3lsT@m8fA{_%|n@Eym=^d++ZHc90knd zx*&I=phm}u8?8o}<0kV^=9pw2${dr;Lz&|w^SCa^y+u%?%yFyLD056P4`q(2=Aq0n z%{-JjPBo9~g4`W~8fA`Gtwx#S74uN$c-cIZIkuaJGRHRarb|=A3_*RH)wf$cQ%c=I zP~Ra^#hq5qvidG5_1%KHz18UM&9j>szJT0PI|`BLf!1a(`h z7f7>-g@XD)s~@s@k(9cPpyt}9&x-{;UnJ=HLP4Jw2>N6Uy)uSg8AE@Jp+CmZD`V)3 zG4#k7dSncJF^0YvGe?V#?TXE=z}rz!I*M_G35ngDhS3@ z6pX1P9uX_VN~>2{{iu}sF+sh{>L;X+izfy3YOB{+{gjmYH9`H1cv?Jb^>bExDfL=G z{kqlbrR&7=f_j718?D|XrQR*5w}{Q+1*>1Q`Xwp#RzdxS)r+O{yhPCRQbC`~1br?S z^voFgV+_4AhF%#%e~h6&#?T{U=#eq>#TYLbLywH1AI8udW9ZEb`e6+HFovEOLr;vM z55~|3W9W%7+^dY?9%jrdf-$cO#_X`aZSb2vOIupMeYK#;f9fBuYVtq%f5pv)}ajhJ zK3>swc)hHvb4G<-S$mCl4W>O?OzV-ah6xSE=N_ z6|?pn@2X~^i=2JkjJLzB0ExLA0HG zy^jCJZRNxzC`sYAAfkk02j2k3;qrl#!dB2U#K$7H`(axOeu0(i9I8HiGrzXXWf6k z3wo{#{<$uU+vJko{*=Ic4W#BphdE-Lh(3p(Eg|2`MSO|BXZbB^uL1-3uu+5Vhu`?JwUef)91eXZZuqJCYQvO61o(eGl{>GiWec<{EI zuG_uY^G{!0%=H@l<5w57_{jCEFzb|FBfoV6_Z50?a_3F1%^5rHZv1Qm*Xivmiq&hi z&~>}-x%D4kn#1)P3i3&^n?#LJ<;0@{HSoPy(K?zZ5kcS{cxo#T&H)N4cePl z)OEZ6sTD;JEp)xo(>i{V=;Qj8eEPmqpZdfNd@tMOl|El@<1at@@yY!;^XT^{ll>0b za!K9>ub<+EEMB#L$Ce3h*p%@%ZGQV!A8GI>`e@vstz$l(`t=kyWLcx!v&(dG!=?^= z@8Fro`ACD`)rWCY-JtDPtsgP@J~w3f+!f!H%BAspcIUdbtHuWner_Me-JkjKncbhv z_3Q1r@Zmc7&f9h4?K<+|y2|D4I`ej2_;8){@pj#KyN-Ogt{(Juo%wKGJnzGG@~OA$ z#@lt|?Yi=Io!R*HM?XGXr+>0P3lHy4xDSrH-{|vJhmAdQUiM^B-!1p3Plw$2p6_+w z)*`dVmhk;&@BBHf^taLTw{JMQ>u33qGWzWT*=LYYe zdHtTRvij1upY!ww{s&+4;X7M?zvkEYb4T=fV`!VDMG_U}_+c$pkEywSzYF6Jx-kAH z`+RmN{u~y@&vnuGus<)`IDMTB`{QhX;`iD9xaf0S=(n)5ex3E}tl#*3zaOvlAsn9! z;}68+?fAs={C;?J-EiB^ScA9E^tjAY0~)N6=Pk~`DtICHK0j z;2*ScxB4*t79ad?d>G#>*2s^>g?W1Tus?Zx*q@JM4Srr9#`W~phqu0b=(lsM!M`jT zm&xOEq4j;Fmle!yM8XS>!!15 zm-xDh#{cK*Cptd);^P?3^R0Fq4#w9_bRC83Av!)&tnbPH(Dl?k#u2*MWw$-l92_KYRIqTV2BR}J~_FD@4+F|d0*XjOAKeldO&9u4; z{sGHd6X!(Zy1qU4?v9gAbCCwWdOR+g|LDg1c)z20qr6BXKlHK2{l1UG`xoYo@}hAA zyp8wvesB91Yw(}4aRcIUc6_|`5v%Qwje9E2iN@J+@sU~|HqK{ye6)Yjyis1Hkspl< z{T}Yaj!Qg$#&N^LckH6?Yq}$ByXbSMy&i36IsE&LxmD*8sp2@1SL72}MK+ONoG7x3 z<3#~cP@Eugh(e;U$SHD(BBH1`Nt`UoigKcuI7O5f{5-d~C?P6}O5#*eQdAaI#A%|G zs4A+7(xQy0E^3IHqLyeR_^#;bqPA!(oH#?&5luu>QCHLx%|x20FB*vEqJ=n9G!!jG zE74lC5#2=((N?q*Jw-3kUhpZrx9B4}icX@h=qEahE~38}Ai9cfVxUMDgT!EQzPLaP z5ktj=;vz9j3>O!RQQ~ZIjuOWO1z+CvFzEi0j1l;#M(5j2AbEsbZR#E@p^(#eL#7al5!*%n>ui9b&GSC+-xp z#C-99xJ%qE7Knx79x+=yC>|1v#KYn-@wiwlmWU_BlVYh@CRU3zV!2o$o)S-sN5o3; zjCfY85|4`KgcobYIR8W ziG$)>@msjHlm88a?^P#6iu@nMkE&C}ajJh3Kda6vvZ?-;_(gSgalAM|uCyA3qUXf1}6Q_v$;zUtgln@0(L2;@mDXNOoL^V-LR2QX14N*qa6lFy% zQBIsL%8S~ff;dA|6m>);QCCzJ^+XlXT+|mWL<7-MoGDs~hN87-B-)6xL|f5Vv=dIW z7fnP5(NuI4%|s`WCi;ucVu0u(28ymCU33$JM0YV*^bkWtPcc;V62nAqFeZc zUOXY57f*@}Vzt;P)`(5wDY02REp~}5;x+MtcwM|Gc8iz98)B zUJ>tzSH-(xhuAB2ihbfe@s)UAd@Vi@--r*ze({kwAU+lc#V6uh@u~Pud?vmZpNk*F z7ve|prGE7-tMn&nHtEmO?9zWpkC*-;Jwf`bG>7z%G^g}8X)cjPaQij)^G2 z1dPOFjKCz6V-BWZHm2f8RA4SDaWtmkD9pfoOvgM_VIgK>0cPP?9D!qSJkG>%I0K7N zjT3Md7ULY8h_i7r&c#Vsf>UumPQiIN9T(y>T!4#k6_(;kT!L$GF|NjCxDJ=%T3n7B zungB@1y*4>Zp2Eg!4+7IwYV2I;U3(A^|%@L;Wj*gTX8?`zy{oo2eA$h<4!z;yYVRQ z!XtPb&*L#XhfR198}S03#LIXBFCm3j@f2Rc)7XN|cn#0u4LpO_@g}z7=lBKQ!rS-| zAK@Lmi(lebcn{m~F@BAo;eGrDpWp**$EWx$euvNSSNsjX$LIJOf5#X265rq-_yhik zf8t;G3V*`C@h$$0zu>#{qki37$~tI-#;A*3&;+}q9`eu>&5)1!*aLf_0d_@mw7_m? zh?Z!Dz0ex%(E)p79~9s~?2G+y5Dv!vNZ=53L>sh4Cmf0c&<>r^1&5(4`k@Hj&>j7e zL=W`D01QMg^u{0z#^LCLAt**)6k;fbp#;M*0TVF-BQXh+QHoKRf^v+;7*t>?#-a?< zP>FFEkLj3^!+*9h<{3s(6TLgu7H3UG@0F~l%sVCPDvH^b9cP_o-Zfco(R(Iq zEsEJ@tY(}w7`=zG{xa{Ntivc~TXvlFnR!QLoks7ctkEcDo3Wa4)^7Cv%DTQaW8pJ>rK!Q#;m(mUEjvSzA(?e?vZ9N z&zN=Bs(p-kF6Y!d`&xId%rj=)@w1OH_k?@M^LY&Si+jwNZLU?55~>)815J2 z#%yz)+Q*n}Jz?Bj&)W;ejoH>4_AzE#0`4Qv>2cgM?lEHtI491dG27e|37G3~+%Lux zu+2HKk1^X^E5^%h4;Pu87#^JLxmu}?>Ig7fNJTR#-et1;`2t7C1fGrGV$`&xHznP<$p<8B{g z-oN*4o_(!5F6J4t?zq{W5k%Z?Rf}t1$691zje?aj#WL>h1X*|54*zaZv*ew9CPhkALf{Aya5`) zW3+?o%p7y=>o}WZu5riFIrA9yvF~p1e4e*48loeNH-Y2Q3HGtCIp*5eF)+tmc;%(2aD=mv9a^H@F679PVq`|baDJSZ_At-B3E1W_4uCnfc`WD0 zdGQztcs_G%^E!IKoCG|UbLkv<4D;-l3)hO*;Bh_oLGT(quW`rA>$S~w;@p~Jn|+;I zbL?Z>@pjB?b4@t6=GbOm=hhti7Cl>=9p{TadPZDhGXOXI>y#r^VW@dUw*cp7uUnR`q{c` z)Vk-jpVw>M(!SoWbxZTSXX{>zdEp)`&2y}*Tbg%JTz7mOQ_pAJy%#L{cb`}{);#`f z-L-C;v9R{7yRL2X{PuHRty`M!yjiz2FU+OqH?JqGTbkD^u3OrtcO3JaVUJl)Ane_+ z*KBVSw}pM=XY0tR23i0h7# z^#Zv54vg!LsdcZ-K3;q1W3Kg(>pHRSy6A+CaBo?6+^lzj>(D)9-SM;T{F~=ohd#Dj z54m2q_3rWCZhe=smHicdvzaQcjsib5o=oTGCd*P$zhrqbUG8jr#z*OXv(`jqDqZF# zcQT%e^HcG4q~hz^>3A52aQ^sj*fQ4>xBuwA)8qNWZ?wYogx^Yqzl96&aGr3TVf?~* zcN~wLa37i9ZQg);osDE!7KE;GLo%KSDc`VA6)BNhMlDfAEJ|NC*`-v%V(zuC_& zGrv_x-!H#Wh<@`G%9*LSe_Vc>X%D{v2xa)2`yKmNWgVZt*^iGO@{7*+@XodcLpFW+ z^rF_=J8w*;KNqKp_HDHOs)p}x`f_cmV!?a)$@J&#sW?B{9Q+gVlck$?KV;$OmC1?K zQB$&V|C2r&wsdQ<^s7?}Hhi=4hj{*=`R$kVsx`lA{oq|c*w`{TCw-h$_*}KAvgC`R zS;;wPMZYObh0mcqD?b&VKNVk3D!#r{Wpw@?Kg6T!TTyF%*cahA@qTzB@9|X?Z?s(y z-7Nov>y7^(%E{bs*gx@cvg7^gude=h@#|lQDz{_Z#&eGs0-`cf}MKNMjQ2BSZc7=mI9z(5SeFpNeC#$Y(cVg$-C65~*c z@fd{)Ou$r3#57DoB_?A!reFrjF&8s&B&u)}j=<5Fg?X5b`Iv(f;QtGYun;HW7%axI zI0?t$WE_vPa0;q%D$d4fI0vU=3C_T|I1?A+JY0hFaVajqWw;Q_uoRc$BCNnQSdOc) z64&7hT#Ku41FppNScO|~BW}hT+=kV-6>D(^Zo=(YhxND<_u+0lfV*%%?!^Y&g9q^t z9>zv&!XtPTPvA*BhR5*~QrL{A@d94NGk6v+;blCB=kW?&#cSAtx9~Pz#~XMD@8V5t I#e3NHJrUY>l>h($ literal 0 HcmV?d00001 diff --git a/science/linear/integration-test/nwp_gal9/resources/nwp_gal9_configuration.nml b/science/linear/integration-test/nwp_gal9/resources/nwp_gal9_configuration.nml new file mode 100644 index 00000000..c687adfe --- /dev/null +++ b/science/linear/integration-test/nwp_gal9/resources/nwp_gal9_configuration.nml @@ -0,0 +1,412 @@ +&base_mesh +file_prefix='resources/mesh_C12_MG', +geometry='spherical', +prepartitioned=.false., +prime_mesh_name='dynamics', +topology='fully_periodic', +/ +&boundaries +limited_area=.false., +transport_overwrite_freq='final', +/ +&checks +limit_cfl=.false., +/ +§ion_choice +aerosol='none', +boundary_layer='none', +chemistry='none', +cloud='none', +dynamics='gungho', +external_forcing=.false., +iau=.false., +iau_sst=.false., +iau_surf=.false., +methane_oxidation=.false., +orographic_drag='none', +radiation='none', +spectral_gwd='none', +stochastic_physics='none', +surface='none', +/ +&convection +dx_ref=50000.0, +l_cvdiag_ctop_qmax=.false., +qlmin=4.0e-4, +resdep_precipramp=.false., +/ +&cosp +l_cosp=.false., +/ +&damping_layer +dl_base=40000.0, +dl_str=0.05, +dl_type='standard', +/ +&departure_points +horizontal_limit='cap', +horizontal_method='ffsl', +n_dep_pt_iterations=1, +share_stencil_extent=.true., +vertical_limit='exponential', +vertical_method='timeaverage', +vertical_sorting=.false., +/ +&energy_correction +encorr_usage='none', +integral_method='fd', +/ +&extrusion +domain_height=80000.0, +method='um_L70_50t_20s_80km', +number_of_layers=70, +planet_radius=6371229.0, +stretching_height=17507.0, +stretching_method='smooth', +/ +&files +ancil_directory='/data/users/lfricadmin/data/ancils/basic-gal/yak/C12', +checkpoint_stem_name='', +diag_stem_name='diagGungho', +ls_directory='/data/users/lfricadmin/data/tangent-linear/Ticket354', +ls_filename='final_ls', +orography_mean_ancil_path='orography/gmted_ramp2/qrparm.orog', +start_dump_directory='/data/users/lfricadmin/data/tangent-linear/Ticket354', +start_dump_filename='final_pert', +/ +&finite_element +cellshape='quadrilateral', +coord_order=1, +coord_system='native', +element_order_h=0, +element_order_v=0, +rehabilitate=.true., +vorticity_in_w1=.false., +/ +&formulation +dlayer_on=.false., +dry_static_adjust=.true., +eos_method='sampled', +exner_from_eos=.false., +horizontal_physics_predictor=.false., +horizontal_transport_predictor=.false., +init_exner_bt=.true., +l_multigrid=.true., +lagged_orog=.true., +moisture_formulation='traditional', +moisture_in_solver=.true., +p2theta_vert=.true., +rotating=.true., +shallow=.true., +si_momentum_equation=.false., +theta_moist_source=.false., +use_multires_coupling=.false., +use_physics=.true., +use_wavedynamics=.true., +vector_invariant=.false., +/ +&helmholtz_solver +gcrk=8, +method='prec_only', +monitor_convergence=.false., +normalise=.true., +preconditioner='multigrid', +si_pressure_a_tol=1.0e-8, +si_pressure_maximum_iterations=400, +si_pressure_tolerance=1.0e-4, +/ +&iau_addinf_io +/ +&iau_addinf_io +/ +&iau_ainc_io +/ +&iau_ainc_io +/ +&iau_bcorr_io +/ +&iau +/ +&idealised +f_lon_deg=0.0, +perturb_init=.false., +test='none', +/ +&ideal_surface +canopy_height=19.01,16.38,0.79,1.26,1.0, +leaf_area_index=5.0,4.0,1.5,1.5,1.5, +n_snow_layers=11*0, +snow_depth=11*0.0, +snow_layer_ice_mass=27*0.0, +snow_layer_temp=27*273.0, +snow_layer_thickness=27*0.0, +soil_moisture=15.86,98.861,274.35,862.27, +soil_temperature=284.508,286.537,289.512,293.066, +surf_tile_fracs=9*0.0,1.0,0.0, +surf_tile_temps=9*295.0,300.0,265.0, +tile_snow_mass=11*0.0, +/ +&initialization +ancil_option='none', +coarse_aerosol_ancil=.false., +coarse_orography_ancil=.false., +coarse_ozone_ancil=.false., +init_option='fd_start_dump', +lbc_option='none', +ls_option='file', +model_eos_height=100, +n_orog_smooth=0, +read_w2h_wind=.true., +sea_ice_source='ancillary', +snow_source='start_dump', +w0_orography_mapping=.false., +zero_w2v_wind=.false., +/ +&initial_density +density_background=0.1, +density_max=2.0, +r1=0.4, +r2=0.4, +x1=0.4, +x2=-0.4, +y1=0.0, +y2=0.0, +z1=0.0, +z2=0.0, +/ +&initial_pressure +method='balanced', +surface_pressure=1000.0e2, +/ +&initial_temperature +bvf_square=0.0001, +pert_centre=60.0, +pert_width_scaling=1.0, +perturb='none', +theta_surf=300.0, +/ +&initial_vapour +/ +&initial_wind +nl_constant=0.0, +profile='constant_uv', +sbr_angle_lat=0.0, +sbr_angle_lon=0.0, +smp_init_wind=.true., +u0=2.0, +v0=0.0, +wind_time_period=0.0, +/ +&io +checkpoint_read=.false., +checkpoint_write=.false., +counter_output_suffix='counter.txt', +diag_active_files='lfric_diag', +diag_always_on_sampling=.false., +diagnostic_frequency=8, +file_convention='UGRID', +multifile_io=.false., +nodal_output_on_w3=.false., +subroutine_counters=.false., +subroutine_timers=.true., +timer_output_path='timer.txt', +use_xios_io=.true., +write_conservation_diag=.false., +write_diag=.false., +write_dump=.false., +write_fluxes=.false., +write_minmax_tseries=.false., +/ +&linear +fixed_ls=.false., +l_stabilise_bl=.false., +ls_read_w2h=.false., +max_bl_stabilisation=0.75, +n_bl_levels_to_stabilise=15, +pert_option='file', +/ +&logging +log_to_rank_zero_only=.false., +run_log_level='debug', +/ +&mixed_solver +eliminate_variables='discrete', +fail_on_non_converged=.true., +gcrk=4, +guess_np1=.false., +mixed_solver_a_tol=1.0e-3, +monitor_convergence=.true., +normalise=.true., +reference_reset_time=1800.0, +si_maximum_iterations=10, +si_method='block_gcr', +si_preconditioner='pressure', +si_tolerance=1.0e-1, +split_w=.true., +/ +&mixing +leonard_term=.false., +smagorinsky=.false., +viscosity=.false., +viscosity_mu=0.0, +/ +&multigrid +chain_mesh_tags='dynamics','multigrid_l1','multigrid_l2', +multigrid_chain_nitems=3, +n_coarsesmooth=4, +n_postsmooth=2, +n_presmooth=2, +smooth_relaxation=0.8, +/ +&esm_couple +l_esm_couple_test=.false., +/ +&orography +orog_init_option='ancil', +/ +&partitioning +generate_inner_halos=.false., +panel_decomposition='auto', +panel_xproc=6, +panel_yproc=1, +partitioner='cubedsphere', +/ +&physics +configure_segments=.false., +limit_drag_incs=.false., +sample_physics_scalars=.true., +sample_physics_winds=.true., +sample_physics_winds_correction=.false., +/ +&planet +cp=1005.0, +gravity=9.80665, +omega=7.292116E-5, +p_zero=100000.0, +rd=287.05, +scaling_factor=1.0, +/ +&radiative_gases +cfc113_rad_opt='off', +cfc11_mix_ratio=1.110e-09, +cfc11_rad_opt='constant', +cfc12_mix_ratio=2.187e-09, +cfc12_rad_opt='constant', +ch4_mix_ratio=1.006e-06, +ch4_rad_opt='constant', +co2_mix_ratio=6.002e-04, +co2_rad_opt='constant', +co_rad_opt='off', +cs_rad_opt='off', +h2_rad_opt='off', +h2o_rad_opt='prognostic', +hcfc22_rad_opt='off', +hcn_rad_opt='off', +he_rad_opt='off', +hfc134a_rad_opt='off', +k_rad_opt='off', +l_cts_fcg_rates=.false., +li_rad_opt='off', +n2_rad_opt='off', +n2o_mix_ratio=4.945e-07, +n2o_rad_opt='constant', +na_rad_opt='off', +nh3_rad_opt='off', +o2_mix_ratio=0.2314, +o2_rad_opt='constant', +o3_rad_opt='ancil', +rb_rad_opt='off', +so2_rad_opt='off', +tio_rad_opt='off', +vo_rad_opt='off', +/ +&solver +gcrk=18, +maximum_iterations=7, +method='chebyshev', +monitor_convergence=.false., +preconditioner='diagonal', +tolerance=1.0e-6, +/ +&specified_surface +/ +&time +calendar='timestep', +calendar_origin='2016-01-01 15:00:00', +calendar_start='2016-01-01 15:00:00', +calendar_type='gregorian', +timestep_end='3', +timestep_start='1', +/ +×tepping +alpha=0.55, +dt=1800, +inner_iterations=2, +method='semi_implicit', +outer_iterations=2, +runge_kutta_method='forward_euler', +spinup_alpha=.false.tau, +tau_r=1.0, +tau_t=1.0, +tau_u=0.55, +/ +&transport +adjust_theta=.false., +adjust_vhv_wind=.false., +ageofair_reset_level=10, +broken_w2_projection=.false., +calculate_detj='upwind', +cap_density_predictor=0.5, +cfl_mol_1d_stab=1.0, +cfl_mol_2d_stab=1.0, +cfl_mol_3d_stab=1.0, +cheap_update=.false., +consistent_metric=.false., +dep_pt_stencil_extent=3, +dry_field_name='density', +enforce_min_value=5*.false., +equation_form=1,2,2,2,2, +ffsl_inner_order=0, +ffsl_outer_order=1, +ffsl_splitting=5*1, +ffsl_unity_3d=.false., +ffsl_vertical_order=2,2,1,2,2, +field_names='density','potential_temperature','wind','moisture', +'con_tracer', +fv_horizontal_order=2, +fv_vertical_order=2, +horizontal_method=5*1, +horizontal_monotone=5*1, +log_space=5*.false., +max_vert_cfl_calc='dep_point', +min_val_abs_tol=-1.0e-12, +min_val_max_iterations=10, +min_val_method='iterative', +min_value=0.0,0.0,-99999999.0,0.0,0.0, +oned_reconstruction=.false., +operators='fv', +panel_edge_high_order=.true., +panel_edge_treatment='none', +profile_size=5, +reversible=5*.false., +runge_kutta_method='ssp3', +scheme=5*1, +si_outer_transport='none', +slice_order='parabola', +special_edges_monotone=5*1, +splitting=5*1, +substep_transport='off', +theta_dispersion_correction=.false., +theta_variable='dry', +transport_ageofair=.false., +use_density_predictor=.false., +vertical_method=5*1, +vertical_monotone=5*1, +vertical_monotone_order=5*3, +vertical_sl_order='cubic', +wind_mono_top=.false., +/ +&validity_test +number_gamma_values=2, +update_ls_frequency=1, +/ diff --git a/science/linear/integration-test/runge_kutta/resources/runge_kutta_configuration.nml b/science/linear/integration-test/runge_kutta/resources/runge_kutta_configuration.nml index 9298f96f..c396ffe9 100644 --- a/science/linear/integration-test/runge_kutta/resources/runge_kutta_configuration.nml +++ b/science/linear/integration-test/runge_kutta/resources/runge_kutta_configuration.nml @@ -131,7 +131,7 @@ subroutine_counters=.false., subroutine_timers=.true., use_xios_io=.false., write_conservation_diag=.false., -write_diag=.true., +write_diag=.false., write_dump=.false., write_fluxes=.false., write_minmax_tseries=.false., @@ -194,7 +194,7 @@ tolerance=1.0e-6, / &time calendar='timestep', -timestep_end='1', +timestep_end='6', timestep_start='1', calendar_type='gregorian', calendar_start='2016-01-01 15:00:00', diff --git a/science/linear/integration-test/runge_kutta/runge_kutta.f90 b/science/linear/integration-test/runge_kutta/runge_kutta.f90 index f460201a..301f0b7f 100644 --- a/science/linear/integration-test/runge_kutta/runge_kutta.f90 +++ b/science/linear/integration-test/runge_kutta/runge_kutta.f90 @@ -10,21 +10,19 @@ !! corresponding nonlinear code. program runge_kutta - use configuration_mod, only: read_configuration, final_configuration use driver_collections_mod, only: init_collections, final_collections use driver_time_mod, only: init_time, final_time + use driver_comm_mod, only: init_comm, final_comm + use driver_log_mod, only: init_logger, final_logger + use driver_config_mod, only: init_config, final_config use driver_modeldb_mod, only: modeldb_type - use halo_comms_mod, only: initialise_halo_comms, & - finalise_halo_comms - use lfric_mpi_mod, only: create_comm, destroy_comm, global_mpi, & - lfric_comm_type - use log_mod, only: initialise_logging, finalise_logging, & - log_event, & + use lfric_mpi_mod, only: global_mpi + use gungho_mod, only: gungho_required_namelists + use log_mod, only: log_event, & LOG_LEVEL_ERROR, & LOG_LEVEL_INFO - use tl_test_driver_mod, only: initialise, & - finalise, & - run_timesteps, & + use linear_driver_mod, only: initialise, finalise + use tl_test_driver_mod, only: run_timesteps_random, & run_kinetic_energy_gradient, & run_advect_density_field, & run_advect_theta_field, & @@ -39,7 +37,6 @@ program runge_kutta ! Model run working data set type(modeldb_type) :: modeldb character(*), parameter :: application_name = "runge_kutta" - character(:), allocatable :: filename ! Variables used for parsing command line arguments @@ -47,8 +44,6 @@ program runge_kutta character(len=0) :: dummy character(len=:), allocatable :: program_name, test_flag - type(lfric_comm_type) :: communicator - ! Flags which determine the tests that will be carried out logical :: do_test_timesteps = .false. logical :: do_test_kinetic_energy_gradient = .false. @@ -65,11 +60,8 @@ program runge_kutta modeldb%mpi => global_mpi - call create_comm( communicator ) - call modeldb%mpi%initialise( communicator ) - call initialise_logging( communicator%get_comm_mpi_val(), & - "linear_integration-runge_kutta-test" ) - call initialise_halo_comms( communicator ) + call modeldb%configuration%initialise( application_name, table_len=10 ) + call modeldb%values%initialise('values', 5) call log_event( 'TL testing running ...', LOG_LEVEL_INFO ) @@ -86,7 +78,7 @@ program runge_kutta call modeldb%fields%add_empty_field_collection("fd_fields", & table_len = 100) - call modeldb%io_contexts%initialise(program_name, 100) + call modeldb%io_contexts%initialise(application_name, 100) ! Parse command line parameters call get_command_argument( 0, dummy, length, status ) @@ -145,16 +137,16 @@ program runge_kutta call log_event( "Unknown test", LOG_LEVEL_ERROR ) end select - call modeldb%configuration%initialise( program_name, table_len=10 ) - call read_configuration( filename, modeldb%configuration ) - deallocate( filename ) - + call init_comm( application_name, modeldb ) + call init_config( filename, gungho_required_namelists, & + modeldb%configuration ) + call init_logger( modeldb%mpi%get_comm(), application_name ) call init_collections() call init_time( modeldb ) - call initialise( application_name, modeldb, modeldb%calendar ) + call initialise( application_name, modeldb ) if (do_test_timesteps) then - call run_timesteps(modeldb) + call run_timesteps_random(modeldb) endif if (do_test_kinetic_energy_gradient) then call run_kinetic_energy_gradient(modeldb) @@ -184,9 +176,8 @@ program runge_kutta call finalise( application_name, modeldb ) call final_time( modeldb ) call final_collections() - call final_configuration() - call finalise_halo_comms() - call finalise_logging() - call destroy_comm() + call final_logger( application_name ) + call final_config() + call final_comm( modeldb ) end program runge_kutta diff --git a/science/linear/integration-test/semi_implicit/resources/semi_implicit_configuration.nml b/science/linear/integration-test/semi_implicit/resources/semi_implicit_configuration.nml index 2b208a45..6be19971 100644 --- a/science/linear/integration-test/semi_implicit/resources/semi_implicit_configuration.nml +++ b/science/linear/integration-test/semi_implicit/resources/semi_implicit_configuration.nml @@ -132,7 +132,7 @@ subroutine_counters=.false., subroutine_timers=.true., use_xios_io=.false., write_conservation_diag=.false., -write_diag=.true., +write_diag=.false., write_dump=.false., write_fluxes=.false., write_minmax_tseries=.false., diff --git a/science/linear/integration-test/semi_implicit/semi_implicit.f90 b/science/linear/integration-test/semi_implicit/semi_implicit.f90 index af625729..57ac5a80 100644 --- a/science/linear/integration-test/semi_implicit/semi_implicit.f90 +++ b/science/linear/integration-test/semi_implicit/semi_implicit.f90 @@ -10,35 +10,25 @@ !! corresponding nonlinear code. program semi_implicit - use configuration_mod, only: read_configuration, final_configuration use driver_collections_mod, only: init_collections, final_collections use driver_time_mod, only: init_time, final_time + use driver_comm_mod, only: init_comm, final_comm + use driver_log_mod, only: init_logger, final_logger + use driver_config_mod, only: init_config, final_config use driver_modeldb_mod, only: modeldb_type - use halo_comms_mod, only: initialise_halo_comms, finalise_halo_comms - use lfric_mpi_mod, only: global_mpi, & - create_comm, destroy_comm, & - lfric_comm_type - use log_mod, only: initialise_logging, finalise_logging, & - log_event, & + use lfric_mpi_mod, only: global_mpi + use gungho_mod, only: gungho_required_namelists + use log_mod, only: log_event, & LOG_LEVEL_ERROR, & LOG_LEVEL_INFO - use namelist_collection_mod, only: namelist_collection_type - use tl_test_driver_mod, only: initialise, & - finalise, & - run_timesteps, & - run_transport_control, & - run_semi_imp_alg, & - run_rhs_sample_eos, & - run_rhs_project_eos, & - run_rhs_alg + use linear_driver_mod, only: initialise, finalise + use tl_test_driver_mod, only: run_timesteps_random implicit none ! Model run working data set type(modeldb_type) :: modeldb - character(*), parameter :: application_name = 'semi_implicit' - character(:), allocatable :: filename ! Variables used for parsing command line arguments @@ -46,26 +36,16 @@ program semi_implicit character(len=0) :: dummy character(len=:), allocatable :: program_name, test_flag - type(lfric_comm_type) :: communicator - ! Flags which determine the tests that will be carried out logical :: do_test_timesteps = .false. - logical :: do_test_transport_control = .false. - logical :: do_test_semi_imp_alg = .false. - logical :: do_test_rhs_alg = .false. - logical :: do_test_rhs_project_eos = .false. - logical :: do_test_rhs_sample_eos = .false. ! Usage message to print character(len=256) :: usage_message modeldb%mpi => global_mpi - call create_comm( communicator ) - call modeldb%mpi%initialise( communicator ) - call initialise_logging( communicator%get_comm_mpi_val(), & - "linear_interface-semi_implicit-test" ) - call initialise_halo_comms( communicator ) + call modeldb%configuration%initialise( application_name, table_len=10 ) + call modeldb%values%initialise('values', 5) call log_event( 'TL testing running ...', LOG_LEVEL_INFO ) @@ -82,7 +62,8 @@ program semi_implicit call modeldb%fields%add_empty_field_collection("fd_fields", & table_len = 100) - call modeldb%io_contexts%initialise(program_name, 100) + + call modeldb%io_contexts%initialise(application_name, 100) ! Parse command line parameters call get_command_argument( 0, dummy, length, status ) @@ -99,8 +80,6 @@ program semi_implicit " transport_control, " // & " semi_imp_alg, " // & " rhs_alg, " // & - " rhs_project_eos, " // & - " rhs_sample_eos, " // & " } " call log_event( trim(usage_message), LOG_LEVEL_ERROR ) end if @@ -118,53 +97,27 @@ program semi_implicit select case (trim(test_flag)) case ("test_timesteps") do_test_timesteps = .true. - case ("test_transport_control") - do_test_transport_control = .true. - case ("test_semi_imp_alg") - do_test_semi_imp_alg = .true. - case ("test_rhs_alg") - do_test_rhs_alg = .true. - case ("test_rhs_project_eos") - do_test_rhs_project_eos = .true. - case ("test_rhs_sample_eos") - do_test_rhs_sample_eos = .true. case default call log_event( "Unknown test", LOG_LEVEL_ERROR ) end select - call modeldb%configuration%initialise( program_name, table_len=10 ) - call read_configuration( filename, modeldb%configuration ) - deallocate( filename ) - + call init_comm( application_name, modeldb ) + call init_config( filename, gungho_required_namelists, & + modeldb%configuration ) + call init_logger( modeldb%mpi%get_comm(), application_name ) call init_collections() call init_time( modeldb ) - call initialise( application_name, modeldb, modeldb%calendar ) + call initialise( application_name, modeldb ) if (do_test_timesteps) then - call run_timesteps(modeldb) - endif - if (do_test_transport_control) then - call run_transport_control(modeldb) - endif - if (do_test_rhs_alg) then - call run_rhs_alg(modeldb) - endif - if (do_test_rhs_project_eos) then - call run_rhs_project_eos(modeldb) - endif - if (do_test_rhs_sample_eos) then - call run_rhs_sample_eos(modeldb) - endif - if (do_test_semi_imp_alg) then - call run_semi_imp_alg(modeldb) + call run_timesteps_random(modeldb) endif call finalise( application_name, modeldb ) call final_time( modeldb ) call final_collections() - call final_configuration() - call finalise_halo_comms() - call finalise_logging() - call destroy_comm() + call final_logger( application_name ) + call final_config() + call final_comm( modeldb ) end program semi_implicit diff --git a/science/linear/integration-test/semi_implicit/semi_implicit.py b/science/linear/integration-test/semi_implicit/semi_implicit.py index af210aff..534eb989 100755 --- a/science/linear/integration-test/semi_implicit/semi_implicit.py +++ b/science/linear/integration-test/semi_implicit/semi_implicit.py @@ -63,43 +63,6 @@ def test_passed(self, out): success = True return success - -class tl_test_semi_imp_alg(TLTest): - ''' - Test the semi implicit timestep - ''' - def __init__(self): - flag = "semi_imp_alg" - super(tl_test_semi_imp_alg, self).__init__(flag) - - -class tl_test_rhs_alg(TLTest): - ''' - Test the right hand side forcing for the mixed solver - ''' - def __init__(self): - flag = "rhs_alg" - super(tl_test_rhs_alg, self).__init__(flag) - -class tl_test_rhs_sample_eos(TLTest): - def __init__(self): - flag = "rhs_sample_eos" - super(tl_test_rhs_sample_eos, self).__init__(flag) - -class tl_test_rhs_project_eos(TLTest): - def __init__(self): - flag = "rhs_project_eos" - super(tl_test_rhs_project_eos, self).__init__(flag) - -class tl_test_transport_control(TLTest): - ''' - Test the transport - ''' - def __init__(self): - flag = "transport_control" - super(tl_test_transport_control, self).__init__(flag) - - class tl_test_timesteps(TLTest): ''' Test running over multiple timesteps @@ -109,9 +72,4 @@ def __init__(self): super(tl_test_timesteps, self).__init__(flag) if __name__ == '__main__': - TestEngine.run( tl_test_rhs_sample_eos() ) - TestEngine.run( tl_test_rhs_project_eos() ) - TestEngine.run(tl_test_transport_control()) - TestEngine.run( tl_test_semi_imp_alg() ) - TestEngine.run( tl_test_rhs_alg() ) - TestEngine.run(tl_test_timesteps()) + TestEngine.run( tl_test_timesteps() ) diff --git a/science/linear/integration-test/tl_test/tl_test_advect_density_field_mod.x90 b/science/linear/integration-test/tl_test/tl_test_advect_density_field_mod.x90 index c7ee7136..b704f347 100644 --- a/science/linear/integration-test/tl_test/tl_test_advect_density_field_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_advect_density_field_mod.x90 @@ -185,7 +185,6 @@ module tl_test_advect_density_field_mod call invoke( X_minus_Y( diff, n2_advection_inc, n1_advection_inc ), & inc_X_minus_Y( diff, p_advection_inc ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) write( log_scratch_space, '(A, E32.12, A, E32.12)' ) & 'gamma = ', gamma_u, ' norm = ', norm_diff diff --git a/science/linear/integration-test/tl_test/tl_test_advect_theta_field_mod.x90 b/science/linear/integration-test/tl_test/tl_test_advect_theta_field_mod.x90 index f712ba30..4630bdf1 100644 --- a/science/linear/integration-test/tl_test/tl_test_advect_theta_field_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_advect_theta_field_mod.x90 @@ -188,8 +188,6 @@ module tl_test_advect_theta_field_mod inc_X_minus_Y( diff, p_advection_inc ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) - write( log_scratch_space, '(A, E32.12, A, E32.12)' ) & 'gamma = ' , gamma_u , ' norm = ' , norm_diff diff --git a/science/linear/integration-test/tl_test/tl_test_convergence_rate_check.f90 b/science/linear/integration-test/tl_test/tl_test_convergence_rate_check.f90 index a13655df..20dee410 100644 --- a/science/linear/integration-test/tl_test/tl_test_convergence_rate_check.f90 +++ b/science/linear/integration-test/tl_test/tl_test_convergence_rate_check.f90 @@ -6,7 +6,7 @@ !>@brief Test the convergence rate (Taylor remainder convergence). module tl_test_convergence_rate_check - use constants_mod, only: r_def, str_def + use constants_mod, only: r_def, str_def, i_def use log_mod, only: log_event, & log_scratch_space, & LOG_LEVEL_INFO @@ -14,17 +14,119 @@ module tl_test_convergence_rate_check implicit none private - public convergence_rate_check + public convergence_pass_string, & + array_convergence_rate_check, & + convergence_rate_check + +contains + + !> @brief Test the convergence rate. + !> @details If the convergence rate is not close to 4, within the + !! specified tolerance, set the pass string to FAIL. + !> @param[in] name Variable or test name + !> @param[in] norm Current norm + !> @param[in] norm_prev Previous norm + !> @param[inout] pass_str Pass string (either PASS or FAIL) + !> @param[in] tol Tolerance + subroutine convergence_pass_string( name, norm, norm_prev, pass_str, tol ) + implicit none + + real(r_def), intent(in) :: norm, norm_prev + character(len=4), intent(inout) :: pass_str + real(r_def), intent(in) :: tol + character(str_def), intent(in) :: name + + real(r_def) :: conv_rate + + conv_rate = sqrt(norm_prev / norm) + + if ( abs( conv_rate - 4.0_r_def ) >= tol ) then + pass_str = "FAIL" + else + pass_str = "PASS" + end if + + write( log_scratch_space, '(A, A, A, E16.8, A, E16.8)') & + name, pass_str, " Convergence rate: ", conv_rate, " Tolerance: ", tol + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + end subroutine convergence_pass_string + + !> @brief Test the convergence rate for individual variables and the sum. + !> @details Calculate the convergence rate based on the norms + !> at two different iterations, and compare with the + !> expected value. Print out either PASS or FAIL, which + !> will then be used by the top level integration test. + !> @param[in] array_norm Norm at second iteration + !> @param[in] array_norm_prev Norm at first iteration + !> @param[in] array_names Name of the variable being tested + !> @param[in] n_variables Array lengths (number of variables) + !> @param[in] label Test name + !> @param[in] tol Tolerance value + subroutine array_convergence_rate_check( array_norm, array_norm_prev, array_names, n_variables, label, tol, indiv_tol) + integer(i_def), intent(in) :: n_variables + real(r_def), intent(in) :: array_norm(n_variables) + real(r_def), intent(in) :: array_norm_prev(n_variables) + character(str_def), intent(in) :: array_names(n_variables) + character(str_def), intent(in) :: label + real(r_def), optional, intent(in) :: tol + real(r_def), optional, intent(in) :: indiv_tol + real(r_def) :: tolerance, individual_tolerance + character(len=4) :: pass_str_arr(n_variables) + character(len=4) :: sum_pass_str, pass_str + character(str_def), parameter :: sum_name = "sum" + integer(i_def) :: i + real(r_def) :: sum, sum_prev + + call log_event( "Checking convergence rate", LOG_LEVEL_INFO ) + + if ( present(tol) ) then + tolerance = tol + else + tolerance = 1.0E-8_r_def + end if + + if ( present(indiv_tol) ) then + individual_tolerance = indiv_tol + else + individual_tolerance = 1.0E-8_r_def + end if + + sum = 0.0_r_def + sum_prev = 0.0_r_def + do i= 1, n_variables + ! Check individual convergence rates + call convergence_pass_string( & + array_names(i), array_norm(i), & + array_norm_prev(i), pass_str_arr(i), & + individual_tolerance ) + + ! Weighted sum + sum = sum + array_norm(i) / array_norm_prev(i) + sum_prev = sum_prev + array_norm_prev(i) / array_norm_prev(i) + end do + + call convergence_pass_string( & + sum_name, sum, sum_prev, sum_pass_str, tolerance) + + pass_str = "PASS" + do i= 1, n_variables + if ( pass_str_arr(i) == "FAIL" ) pass_str = "FAIL" + end do + if ( sum_pass_str == "FAIL" ) pass_str = "FAIL" + + write(log_scratch_space,'(" test",A32," : ",A4)') trim(label), pass_str + call log_event( log_scratch_space, LOG_LEVEL_INFO ) - contains + end subroutine array_convergence_rate_check !> @brief Calculate and test the convergence rate. !> @details Calculate the convergence rate based on the norms !> at two different iterations, and compare with the !> expected value. Print out either PASS or FAIL, which !> will then be used by the top level integration test. - !> @param[in] norm_diff Norm at first iteration - !> @param[in] norm_diff_prev Norm at second iteration + !> @param[in] norm_diff Norm at second iteration + !> @param[in] norm_diff_prev Norm at first iteration !> @param[in] label Name of the code being tested !> @param[in] tol Tolerance value subroutine convergence_rate_check( norm_diff, norm_diff_prev, label, tol ) @@ -35,7 +137,6 @@ subroutine convergence_rate_check( norm_diff, norm_diff_prev, label, tol ) character(str_def), intent(in) :: label real(r_def), optional, intent(in) :: tol real(r_def) :: tolerance - real(r_def) :: conv_rate character(len=4) :: pass_str if ( present(tol) ) then @@ -44,23 +145,8 @@ subroutine convergence_rate_check( norm_diff, norm_diff_prev, label, tol ) tolerance = 1.0E-8_r_def end if - conv_rate = norm_diff_prev/ norm_diff - - - - write( log_scratch_space, '(A)' ) & - "TL Test: " // trim(label) - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - write( log_scratch_space, '(A, E16.8)') & - "Convergence rate: ", conv_rate - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - - if ( abs(conv_rate - 4.0_r_def ) < tolerance ) then - pass_str = "PASS" - else - pass_str = "FAIL" - end if + call convergence_pass_string( & + label, norm_diff, norm_diff_prev, pass_str, tolerance) write(log_scratch_space,'(" test",A32," : ",A4)') trim(label), pass_str call log_event( log_scratch_space, LOG_LEVEL_INFO ) diff --git a/science/linear/integration-test/tl_test/tl_test_driver_mod.f90 b/science/linear/integration-test/tl_test/tl_test_driver_mod.f90 index b93d850d..0d1c7e5e 100644 --- a/science/linear/integration-test/tl_test/tl_test_driver_mod.f90 +++ b/science/linear/integration-test/tl_test/tl_test_driver_mod.f90 @@ -48,13 +48,13 @@ module tl_test_driver_mod use tl_test_rhs_alg_mod, only : test_rhs_alg use tl_test_semi_imp_alg_mod, only : test_semi_imp_alg use tl_test_timesteps_alg_mod, only : test_timesteps + use tl_test_timesteps_random_alg_mod, only : test_timesteps_random implicit none private - public initialise, & - finalise, & - run_timesteps, & + public run_timesteps, & + run_timesteps_random, & run_kinetic_energy_gradient, & run_advect_density_field, & run_advect_theta_field, & @@ -72,109 +72,47 @@ module tl_test_driver_mod type(mesh_type), pointer :: mesh => null() type(mesh_type), pointer :: twod_mesh => null() - type(mesh_type), pointer :: aerosol_mesh => null() - type(mesh_type), pointer :: aerosol_twod_mesh => null() contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !>@brief Sets up the required state in preparation for run. - !>@param [in] program_name An identifier given to the model being run - !>@param [in,out] modeldb The structure that holds model state - !> - subroutine initialise( program_name, modeldb, calendar ) + !>@brief Tests the tangent linear model for multiple timesteps + !>@param [in,out] modeldb The structure that holds model state + subroutine run_timesteps(modeldb) implicit none - character(*), intent(in) :: program_name - type(modeldb_type), intent(inout) :: modeldb - class(calendar_type), intent(in) :: calendar - - type(gungho_time_axes_type) :: model_axes - type(io_value_type) :: temp_corr_io_value - type(io_value_type) :: random_seed_io_value - integer(i_def) :: random_seed_size - real(r_def), allocatable :: real_array(:) - - call modeldb%values%initialise( 'values', 5 ) - - ! Initialise infrastructure and setup constants - ! - call initialise_infrastructure( program_name, modeldb ) - - ! Add a place to store time axes in modeldb - call modeldb%values%add_key_value('model_axes', model_axes) + type(modeldb_type), intent(inout) :: modeldb - ! Get primary and 2D meshes for initialising model data mesh => mesh_collection%get_mesh(prime_mesh_name) twod_mesh => mesh_collection%get_mesh(mesh, TWOD) - ! Assume aerosol mesh is the same as dynamics mesh - aerosol_mesh => mesh - aerosol_twod_mesh => twod_mesh - - ! gungho_init_field() expects these values to exist. The dependency of - ! the linear application tests on this procedure will hopefully be resolved - ! in the future, at which point this initialisation may be removed. - ! - call temp_corr_io_value%init('temperature_correction_rate', [0.0_r_def]) - call modeldb%values%add_key_value( 'temperature_correction_io_value', & - temp_corr_io_value ) - call modeldb%values%add_key_value( 'total_dry_mass', 0.0_r_def ) - call modeldb%values%add_key_value( 'total_energy', 0.0_r_def ) - call modeldb%values%add_key_value( 'total_energy_previous', 0.0_r_def ) - if ( stochastic_physics == stochastic_physics_um ) then - ! Random seed for stochastic physics - call random_seed(size = random_seed_size) - allocate(real_array(random_seed_size)) - real_array(1:random_seed_size) = 0.0_r_def - call random_seed_io_value%init("random_seed", real_array) - call modeldb%values%add_key_value( 'random_seed_io_value', & - random_seed_io_value ) - deallocate(real_array) - end if - - ! Instantiate the fields stored in model_data - call create_model_data( modeldb, & - mesh, & - twod_mesh, & - aerosol_mesh, & - aerosol_twod_mesh ) - - ! Instantiate the linearisation state - call linear_create_ls( modeldb, mesh ) - - ! Initialise the fields stored in the model_data prognostics. This needs - ! to be done before initialise_model. - call initialise_model_data( modeldb, mesh, twod_mesh ) - - ! Model configuration initialisation - call initialise_model( mesh, & - modeldb ) - - ! Initialise the linearisation state - call linear_init_ls( mesh, twod_mesh, modeldb ) - - ! Finalise model - call finalise_model(modeldb) - - end subroutine initialise + call test_timesteps( modeldb, & + mesh, & + twod_mesh, & + modeldb%clock ) + + end subroutine run_timesteps !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !>@brief Tests the tangent linear model for multiple timesteps + !! using prescribed random data for the initial conditions !>@param [in,out] modeldb The structure that holds model state - subroutine run_timesteps(modeldb) + subroutine run_timesteps_random(modeldb) implicit none type(modeldb_type), intent(inout) :: modeldb - call test_timesteps( modeldb, & - mesh, & - twod_mesh, & - modeldb%clock ) + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) - end subroutine run_timesteps + call test_timesteps_random( modeldb, & + mesh, & + twod_mesh, & + modeldb%clock ) + + end subroutine run_timesteps_random !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !>@brief Tests the tangent linear model kinetic energy gradient kernel @@ -185,6 +123,9 @@ subroutine run_kinetic_energy_gradient(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_kinetic_energy_gradient( modeldb, & mesh, & twod_mesh ) @@ -200,6 +141,9 @@ subroutine run_advect_density_field(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_advect_density_field( modeldb, & mesh, & twod_mesh ) @@ -215,6 +159,9 @@ subroutine run_advect_theta_field(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_advect_theta_field( modeldb, & mesh, & twod_mesh ) @@ -230,6 +177,9 @@ subroutine run_vorticity_advection(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_vorticity_advection( modeldb, & mesh, & twod_mesh ) @@ -245,6 +195,9 @@ subroutine run_project_eos_pressure(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_project_eos_pressure( modeldb, & mesh, & twod_mesh ) @@ -260,6 +213,9 @@ subroutine run_sample_eos_pressure(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_sample_eos_pressure( modeldb, & mesh, & twod_mesh ) @@ -275,6 +231,9 @@ subroutine run_hydrostatic(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_hydrostatic( modeldb, & mesh, & twod_mesh ) @@ -290,6 +249,9 @@ subroutine run_pressure_gradient_bd(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_pressure_gradient_bd( modeldb, & mesh, & twod_mesh ) @@ -305,6 +267,8 @@ subroutine run_rk_alg(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + call test_rk_alg( modeldb, & mesh) @@ -319,6 +283,9 @@ subroutine run_transport_control(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_transport_control( modeldb, & mesh, & twod_mesh ) @@ -334,6 +301,9 @@ subroutine run_semi_imp_alg(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_semi_imp_alg( modeldb, & mesh, & twod_mesh ) @@ -349,6 +319,9 @@ subroutine run_rhs_alg(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_rhs_alg( modeldb, & mesh, & twod_mesh ) @@ -364,6 +337,9 @@ subroutine run_rhs_project_eos(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_rhs_project_eos( modeldb, & mesh, & twod_mesh ) @@ -379,31 +355,13 @@ subroutine run_rhs_sample_eos(modeldb) type(modeldb_type), intent(inout) :: modeldb + mesh => mesh_collection%get_mesh(prime_mesh_name) + twod_mesh => mesh_collection%get_mesh(mesh, TWOD) + call test_rhs_sample_eos( modeldb, & mesh, & twod_mesh ) end subroutine run_rhs_sample_eos - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !>@brief Tidies up after a run. - !>@param [in] program_name An identifier given to the model being run - !>@param [in,out] modeldb The structure that holds model state - subroutine finalise( program_name, modeldb ) - - implicit none - - character(*), intent(in) :: program_name - type(modeldb_type), intent(inout) :: modeldb - - call log_event( 'Finalising '//program_name//' ...', LOG_LEVEL_ALWAYS ) - - ! Destroy the fields stored in model_data - call finalise_model_data( modeldb ) - - ! Finalise infrastructure and constants - call finalise_infrastructure(modeldb) - - end subroutine finalise - end module tl_test_driver_mod diff --git a/science/linear/integration-test/tl_test/tl_test_hydrostatic_mod.x90 b/science/linear/integration-test/tl_test/tl_test_hydrostatic_mod.x90 index fbcb7baa..5aa4f198 100644 --- a/science/linear/integration-test/tl_test/tl_test_hydrostatic_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_hydrostatic_mod.x90 @@ -189,8 +189,6 @@ module tl_test_hydrostatic_mod inc_X_minus_Y( diff, p_rhs_u ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) - write( log_scratch_space, '(A, E32.12, A, E32.12)' ) & 'gamma = ', gamma, ' norm = ', norm_diff diff --git a/science/linear/integration-test/tl_test/tl_test_kinetic_energy_gradient_mod.x90 b/science/linear/integration-test/tl_test/tl_test_kinetic_energy_gradient_mod.x90 index 2138662a..f2910b1b 100644 --- a/science/linear/integration-test/tl_test/tl_test_kinetic_energy_gradient_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_kinetic_energy_gradient_mod.x90 @@ -135,7 +135,6 @@ module tl_test_kinetic_energy_gradient_mod call invoke( X_minus_Y( diff, n2_rhs_u, n1_rhs_u ), & inc_X_minus_Y( diff, p_rhs_u ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) write( log_scratch_space, '(A, E32.12, A, E32.12)' ) & 'gamma = ', gamma, ' norm = ', norm_diff diff --git a/science/linear/integration-test/tl_test/tl_test_pressure_grad_bd_mod.x90 b/science/linear/integration-test/tl_test/tl_test_pressure_grad_bd_mod.x90 index 8bb6117a..5a9f0210 100644 --- a/science/linear/integration-test/tl_test/tl_test_pressure_grad_bd_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_pressure_grad_bd_mod.x90 @@ -213,8 +213,6 @@ module tl_test_pressure_grad_bd_mod inc_X_minus_Y( diff, p_rhs_u ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) - write( log_scratch_space, '(A, E32.12, A, E32.12)' ) & 'gamma = ', gamma, ' norm = ', norm_diff diff --git a/science/linear/integration-test/tl_test/tl_test_project_eos_pressure_mod.x90 b/science/linear/integration-test/tl_test/tl_test_project_eos_pressure_mod.x90 index 01ef9b69..9fc26e06 100644 --- a/science/linear/integration-test/tl_test/tl_test_project_eos_pressure_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_project_eos_pressure_mod.x90 @@ -195,7 +195,6 @@ module tl_test_project_eos_pressure_mod call invoke( X_minus_Y( diff, n2_exner, n1_exner ), & inc_X_minus_Y( diff, p_exner ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) write( log_scratch_space, '(A, E32.12, A, E32.12)' ) & 'gamma = ', gamma , ' norm = ', norm_diff diff --git a/science/linear/integration-test/tl_test/tl_test_rhs_alg_mod.x90 b/science/linear/integration-test/tl_test/tl_test_rhs_alg_mod.x90 index 1c483372..794764e8 100644 --- a/science/linear/integration-test/tl_test/tl_test_rhs_alg_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_rhs_alg_mod.x90 @@ -18,7 +18,7 @@ module tl_test_rhs_alg_mod use function_space_mod, only: function_space_type use derived_config_mod, only: bundle_size use moist_dyn_mod, only: num_moist_factors - use sci_assign_field_random_kernel_mod, only: assign_field_random_kernel_type + use mr_indices_mod, only: nummr use field_indices_mod, only: igh_u, igh_t, igh_d, igh_p use sci_field_bundle_builtins_mod, only: clone_bundle, & add_bundle, & @@ -28,8 +28,10 @@ module tl_test_rhs_alg_mod use rhs_alg_mod, only: rhs_alg use tl_rhs_alg_mod, only: tl_rhs_alg use log_mod, only: log_event, LOG_LEVEL_INFO, & - LOG_LEVEL_ERROR - use tl_test_convergence_rate_check, only: convergence_rate_check + LOG_LEVEL_ERROR + use moist_dyn_factors_alg_mod, only: moist_dyn_factors_alg + use tl_moist_dyn_factors_alg_mod, only: tl_moist_dyn_factors_alg + use tl_test_convergence_rate_check, only: array_convergence_rate_check implicit none @@ -56,8 +58,8 @@ module tl_test_rhs_alg_mod character(str_def) :: label = "rhs_alg" - type(field_collection_type ), pointer :: ls_fields - + type( field_collection_type ), pointer :: ls_fields + type( field_collection_type ), pointer :: prognostic_fields => null() type( field_type ) :: state(bundle_size) type( field_type ) :: ls_state(bundle_size) @@ -73,37 +75,62 @@ module tl_test_rhs_alg_mod type(field_type), pointer :: ls_theta => null() type(field_type), pointer :: ls_exner => null() type(field_type), pointer :: ls_moist_dyn(:) => null() + type(field_type), pointer :: ls_mr(:) => null() + type(field_type), pointer :: u => null() + type(field_type), pointer :: rho => null() + type(field_type), pointer :: theta => null() + type(field_type), pointer :: exner => null() + type(field_type), pointer :: moist_dyn(:) => null() + type(field_type), pointer :: mr(:) => null() type(field_type), dimension(num_moist_factors) :: p_moist_dyn type(field_type), dimension(num_moist_factors) :: n_moist_dyn type(field_type), dimension(num_moist_factors) :: r_moist_dyn type(field_collection_type), pointer :: moisture_fields => null() + type(field_array_type), pointer :: mr_array => null() + type(field_array_type), pointer :: moist_dyn_array => null() + type(field_array_type), pointer :: ls_mr_array => null() type(field_array_type), pointer :: ls_moist_dyn_array => null() real(r_def) :: gamma_u, gamma_rho, gamma_exner, gamma_theta, & gamma_moist_dyn - real(r_def) :: norm_u, norm_rho, norm_exner, norm_theta, norm_moist_dyn - real(r_def) :: norm_moist_dyn_tmp - real(r_def) :: norm_diff, norm_diff_prev + real(r_def) :: norm_u, norm_exner + integer(i_def), parameter :: n_variables = 2 + real(r_def) :: array_norm(n_variables), array_norm_prev(n_variables) + character(str_def) :: array_names(n_variables) - real(r_def), parameter :: tol = 1.e-2_r_def + real(r_def), parameter :: tol = 5.e-2_r_def + real(r_def), parameter :: indiv_tol = 5.e-1_r_def integer :: i, n call log_event( "TL Test: " // trim(label), & LOG_LEVEL_INFO ) + prognostic_fields => modeldb%fields%get_field_collection( & + "prognostic_fields") ls_fields => modeldb%fields%get_field_collection("ls_fields") moisture_fields => modeldb%fields%get_field_collection("moisture_fields") + call moisture_fields%get_field("mr", mr_array) + call moisture_fields%get_field("moist_dyn", moist_dyn_array) + mr => mr_array%bundle + moist_dyn => moist_dyn_array%bundle + call moisture_fields%get_field("ls_mr", ls_mr_array) call moisture_fields%get_field("ls_moist_dyn", ls_moist_dyn_array) + ls_mr => ls_mr_array%bundle ls_moist_dyn => ls_moist_dyn_array%bundle - ! Input + ! Linearisation state call ls_fields%get_field('ls_u', ls_u) call ls_fields%get_field('ls_rho', ls_rho) call ls_fields%get_field('ls_theta', ls_theta) call ls_fields%get_field('ls_exner', ls_exner) + ! Perturbation + call prognostic_fields%get_field('u', u) + call prognostic_fields%get_field('rho', rho) + call prognostic_fields%get_field('theta', theta) + call prognostic_fields%get_field('exner', exner) call ls_state(igh_u)%initialise( vector_space = ls_u%get_function_space() ) call ls_state(igh_t)%initialise( vector_space = ls_theta%get_function_space() ) @@ -127,24 +154,27 @@ module tl_test_rhs_alg_mod setval_X(ls_state(igh_d), ls_rho ), & setval_X(ls_state(igh_p), ls_exner) ) - call invoke( assign_field_random_kernel_type( random(igh_u), 1.0_r_def ) , & - assign_field_random_kernel_type( random(igh_d), 1.0_r_def ) , & - assign_field_random_kernel_type( random(igh_t), 1.0_r_def ) , & - assign_field_random_kernel_type( random(igh_p), 1.0_r_def ) ) + call moist_dyn_factors_alg(ls_moist_dyn, ls_mr) - do i = 1, num_moist_factors - call invoke( assign_field_random_kernel_type( r_moist_dyn(i), 1.0_r_def ) ) - enddo + ! Set the 'random data' to be the perturbation + call invoke( name = "copy_fields_to_state", & + setval_X(random(igh_u), u ), & + setval_X(random(igh_p), exner ), & + setval_X(random(igh_t), theta), & + setval_X(random(igh_d), rho ) ) - gamma_u = 1.e2_r_def - gamma_theta = 1.e2_r_def - gamma_rho = 1.e2_r_def - gamma_exner = 1.e1_r_def - gamma_moist_dyn = 1.e-1_r_def + call tl_moist_dyn_factors_alg(r_moist_dyn, mr ) + + gamma_u = 1.0_r_def + gamma_theta = 1.0_r_def + gamma_rho = 1.0_r_def + gamma_exner = 1.0_r_def + gamma_moist_dyn = 1.0_r_def call set_bundle_scalar( 0.0_r_def, n1_rhs, bundle_size ) call rhs_alg( n1_rhs, dt, ls_state, ls_state, ls_moist_dyn, & - .true., .true., .false., modeldb%clock ) + compute_eos=.true., compute_rhs_t_d=.true., & + dlayer_rhs=.true., model_clock=modeldb%clock ) do n=1,2 gamma_u = gamma_u / 2.0_r_def @@ -166,42 +196,45 @@ module tl_test_rhs_alg_mod do i = 1, num_moist_factors call invoke( & a_times_X( p_moist_dyn(i), gamma_moist_dyn, r_moist_dyn(i) ), & - setval_X( n_moist_dyn(i), ls_moist_dyn(i) ), & + setval_X( n_moist_dyn(i), ls_moist_dyn(i) ), & inc_X_plus_Y( n_moist_dyn(i), p_moist_dyn(i) ) ) enddo call rhs_alg( n2_rhs, dt, state, state, n_moist_dyn, & - .true., .true., .false., modeldb%clock ) + compute_eos=.true., compute_rhs_t_d=.true., & + dlayer_rhs=.true., model_clock=modeldb%clock ) call tl_rhs_alg(p_rhs, dt, p_state, p_state, p_moist_dyn, & ls_state, ls_moist_dyn, & - .true., .false., modeldb%clock) + compute_eos=.true., dlayer_rhs=.true., & + model_clock=modeldb%clock) ! diff = n2_rhs - n1_rhs call minus_bundle( n2_rhs, n1_rhs, diff, bundle_size ) call invoke( & inc_X_minus_Y( diff(igh_u), p_rhs(igh_u) ), & - inc_X_minus_Y( diff(igh_t), p_rhs(igh_t) ), & - inc_X_minus_Y( diff(igh_d), p_rhs(igh_d) ), & inc_X_minus_Y( diff(igh_p), p_rhs(igh_p) ) ) call invoke( X_innerproduct_X( norm_u, diff(igh_u) ) , & - X_innerproduct_X( norm_rho, diff(igh_d) ) , & - X_innerproduct_X( norm_theta, diff(igh_t) ) , & X_innerproduct_X( norm_exner, diff(igh_p) ) ) - print*, norm_u, norm_rho, norm_theta, norm_exner - norm_diff = norm_u + norm_rho + norm_theta + norm_exner - norm_diff = sqrt(norm_diff) - print*,'norm', norm_diff + ! The equations for rho and theta are already linear + ! so only evaluate u and exner + array_norm(1) = norm_u + array_norm(2) = norm_exner + array_names(1) = 'norm_u' + array_names(2) = 'norm_exner' if (n == 2) then - call convergence_rate_check( norm_diff, norm_diff_prev, & - label, tol=tol ) - endif + call array_convergence_rate_check( & + array_norm, array_norm_prev, & + array_names, n_variables, label, tol=tol, & + indiv_tol=indiv_tol ) + end if + + array_norm_prev = array_norm - norm_diff_prev = norm_diff enddo end subroutine test_rhs_alg diff --git a/science/linear/integration-test/tl_test/tl_test_rhs_project_eos_mod.x90 b/science/linear/integration-test/tl_test/tl_test_rhs_project_eos_mod.x90 index 9504122f..8b39ecb5 100644 --- a/science/linear/integration-test/tl_test/tl_test_rhs_project_eos_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_rhs_project_eos_mod.x90 @@ -209,7 +209,6 @@ contains call invoke( X_minus_Y( diff, n2_rhs, n1_rhs ), & inc_X_minus_Y( diff, p_rhs ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) if (n == 2) then call convergence_rate_check( norm_diff, norm_diff_prev, & diff --git a/science/linear/integration-test/tl_test/tl_test_rhs_sample_eos_mod.x90 b/science/linear/integration-test/tl_test/tl_test_rhs_sample_eos_mod.x90 index 49d69273..0d1cd26d 100644 --- a/science/linear/integration-test/tl_test/tl_test_rhs_sample_eos_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_rhs_sample_eos_mod.x90 @@ -195,7 +195,6 @@ contains call invoke( X_minus_Y( diff, n2_rhs, n1_rhs ), & inc_X_minus_Y( diff, p_rhs ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) if (n == 2) then call convergence_rate_check( norm_diff, norm_diff_prev, & diff --git a/science/linear/integration-test/tl_test/tl_test_rk_alg_mod.x90 b/science/linear/integration-test/tl_test/tl_test_rk_alg_mod.x90 index 9b99f0fc..d16f4c14 100644 --- a/science/linear/integration-test/tl_test/tl_test_rk_alg_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_rk_alg_mod.x90 @@ -251,6 +251,7 @@ module tl_test_rk_alg_mod ! Evolve the perturbation fields with the linear timestep. This ! gives p_x=L(gamma r_x) + call tl_rk_alg_final() call tl_rk_alg_init(mesh, p_u, p_rho, p_theta, p_exner, & ls_u, ls_rho, ls_theta, ls_exner) @@ -300,7 +301,6 @@ module tl_test_rk_alg_mod call log_event( log_scratch_space, LOG_LEVEL_INFO ) norm_diff = norm_rho + norm_theta + norm_exner + norm_u + norm_moist_dyn - norm_diff = sqrt(norm_diff) if (n == 2) then ! Compare results from first and second iteration diff --git a/science/linear/integration-test/tl_test/tl_test_sample_eos_pressure_mod.x90 b/science/linear/integration-test/tl_test/tl_test_sample_eos_pressure_mod.x90 index b6f56add..1915214c 100644 --- a/science/linear/integration-test/tl_test/tl_test_sample_eos_pressure_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_sample_eos_pressure_mod.x90 @@ -175,7 +175,6 @@ module tl_test_sample_eos_pressure_mod call invoke( X_minus_Y( diff, n2_exner, n1_exner ), & inc_X_minus_Y( diff, p_exner ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) write( log_scratch_space, '(A, E32.12, A, E32.12)' ) & 'gamma = ', gamma , ' norm = ', norm_diff diff --git a/science/linear/integration-test/tl_test/tl_test_semi_imp_alg_mod.x90 b/science/linear/integration-test/tl_test/tl_test_semi_imp_alg_mod.x90 index 2b060eac..6a58e383 100644 --- a/science/linear/integration-test/tl_test/tl_test_semi_imp_alg_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_semi_imp_alg_mod.x90 @@ -10,8 +10,6 @@ !! implicit algorithm, using the Taylor remainder convergence test. module tl_test_semi_imp_alg_mod - use sci_assign_field_random_kernel_mod, & - only: assign_field_random_kernel_type use constants_mod, only: i_def, r_def, str_def use field_array_mod, only: field_array_type use field_mod, only: field_type @@ -39,7 +37,7 @@ module tl_test_semi_imp_alg_mod use timestep_method_mod, only: timestep_method_type use moist_dyn_factors_alg_mod, only: moist_dyn_factors_alg use tl_moist_dyn_factors_alg_mod, only: tl_moist_dyn_factors_alg - use tl_test_convergence_rate_check, only: convergence_rate_check + use tl_test_convergence_rate_check, only: array_convergence_rate_check implicit none @@ -128,9 +126,12 @@ module tl_test_semi_imp_alg_mod real(r_def) :: gamma_mr real(r_def) :: norm_u, norm_rho, norm_exner, norm_theta real(r_def) :: norm_mr, norm_mr_tmp - real(r_def) :: norm_diff, norm_diff_prev + integer(i_def), parameter :: n_variables = 5 + real(r_def) :: array_norm(n_variables), array_norm_prev(n_variables) + character(str_def) :: array_names(n_variables) - real(r_def), parameter :: tol = 5.e-1_r_def + real(r_def), parameter :: tol = 5.e-2_r_def + real(r_def), parameter :: indiv_tol = 5.e-1_r_def real(r_def), parameter :: dtemp_encorr = 0.0_r_def integer :: n, i @@ -173,12 +174,14 @@ module tl_test_semi_imp_alg_mod aerosol_fields => modeldb%fields%get_field_collection("aerosol_fields") stph_fields => modeldb%fields%get_field_collection("stph_fields") lbc_fields => modeldb%fields%get_field_collection("lbc_fields") - ! Input + + ! Linearisation State call ls_fields%get_field('ls_u', ls_u) call ls_fields%get_field('ls_rho', ls_rho) call ls_fields%get_field('ls_theta', ls_theta) call ls_fields%get_field('ls_exner', ls_exner) + ! Perturbation call prognostic_fields%get_field('u', u) call prognostic_fields%get_field('rho', rho) call prognostic_fields%get_field('theta', theta) @@ -228,15 +231,23 @@ module tl_test_semi_imp_alg_mod call clone_bundle(mr, b_mr, nummr) call clone_bundle(mr, diff_mr, nummr) - call invoke( assign_field_random_kernel_type( r_u, 1.0_r_def ) , & - assign_field_random_kernel_type( r_rho, 1.0_r_def ) , & - assign_field_random_kernel_type( r_theta, 1.0_r_def ) , & - assign_field_random_kernel_type( r_exner, 1.0_r_def ) ) + ! Set the 'random data' to be the perturbation + call invoke( name = "copy_fields_to_state", & + setval_X(r_u, u ), & + setval_X(r_exner, exner ), & + setval_X(r_theta, theta), & + setval_X(r_rho, rho ) ) + ! Overwrite the moisture fields with potential temperature do i = 1, nummr - call invoke( assign_field_random_kernel_type( r_mr(i), 1.0_r_def ) ) + call invoke( & + setval_X( r_mr(i), theta ), & + inc_a_times_X( 0.000001, r_mr(i) ), & + setval_X( ls_mr(i), ls_theta ), & + inc_a_times_X( 0.000001, ls_mr(i) )) end do + ! Set the prognostic fields to be the linearisation state call invoke( setval_X( u, ls_u ), & setval_X( theta, ls_theta ), & setval_X( rho, ls_rho ), & @@ -247,7 +258,9 @@ module tl_test_semi_imp_alg_mod end do do i = 1, nummr call invoke( setval_X( mr(i), ls_mr(i) ) ) - end do + end do + + call moist_dyn_factors_alg(moist_dyn, mr) ! Clean up solvers left over from previous tests, as semi-implicit ! method will create new solvers which cause double-allocates @@ -274,15 +287,15 @@ module tl_test_semi_imp_alg_mod call invoke( setval_X( b_mr(i), mr(i) ) ) end do - gamma_u = 1.e4_r_def - gamma_theta = 1.e2_r_def - gamma_rho = 8.0e-1_r_def - gamma_exner = 6.5e-1_r_def + gamma_u = 6.5_r_def + gamma_theta = 6.5_r_def + gamma_rho = 6.5_r_def + gamma_exner = 6.5_r_def if ( moisture_formulation == moisture_formulation_dry ) then gamma_mr = 0.0_r_def else - gamma_mr = 1.e-2_r_def + gamma_mr = 6.5_r_def end if do n=1,2 @@ -328,6 +341,8 @@ module tl_test_semi_imp_alg_mod ! Evolve the perturbation fields with the linear timestep. This ! gives p_x=L(gamma r_x) + call final_si_operators() + call tl_semi_implicit_alg_final() call tl_semi_implicit_alg_init(mesh, p_u, p_rho, p_theta, p_exner, & p_mr, ls_u, ls_rho, ls_theta, ls_exner, & ls_mr, ls_moist_dyn) @@ -365,35 +380,24 @@ module tl_test_semi_imp_alg_mod norm_mr = norm_mr + norm_mr_tmp end do - write( log_scratch_space, & - '(A, E32.12, A, E32.12 )' ) & - ' norm_u = ' , norm_u, & - ' norm_rho = ' , norm_rho - - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - write( log_scratch_space, & - '(A, E32.12, A, E32.12 )' ) & - ' norm_mr = ' , norm_mr, & - ' norm_exner = ' , norm_exner - - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - write( log_scratch_space, & - '(A, E32.12 )' ) & - ' norm_theta = ' , norm_theta - - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - norm_diff = norm_rho + norm_theta + norm_exner + norm_u + norm_mr - norm_diff = sqrt(norm_diff) + array_norm(1) = norm_u + array_norm(2) = norm_rho + array_norm(3) = norm_theta + array_norm(4) = norm_mr + array_norm(5) = norm_exner + array_names(1) = 'norm_u' + array_names(2) = 'norm_rho' + array_names(3) = 'norm_theta' + array_names(4) = 'norm_mr' + array_names(5) = 'norm_exner' if (n == 2) then - call convergence_rate_check( norm_diff, norm_diff_prev, & - label, tol=tol ) + call array_convergence_rate_check( array_norm, array_norm_prev, & + array_names, n_variables, label, tol=tol, indiv_tol=indiv_tol ) end if - norm_diff_prev = norm_diff + array_norm_prev = array_norm + end do end subroutine test_semi_imp_alg diff --git a/science/linear/integration-test/tl_test/tl_test_timesteps_alg_mod.x90 b/science/linear/integration-test/tl_test/tl_test_timesteps_alg_mod.x90 index 65759abf..7e5552de 100644 --- a/science/linear/integration-test/tl_test/tl_test_timesteps_alg_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_timesteps_alg_mod.x90 @@ -44,7 +44,7 @@ module tl_test_timesteps_alg_mod method, & method_rk use time_config_mod, only: timestep_start, timestep_end - use tl_test_convergence_rate_check, only: convergence_rate_check + use tl_test_convergence_rate_check, only: array_convergence_rate_check use validity_test_config_mod, only: number_gamma_values, & update_ls_frequency use moist_dyn_factors_alg_mod, only: moist_dyn_factors_alg @@ -52,7 +52,6 @@ module tl_test_timesteps_alg_mod implicit none - private convergence_pass_string public test_timesteps contains @@ -120,8 +119,12 @@ module tl_test_timesteps_alg_mod real(r_def) :: denom_u, denom_rho, denom_exner, denom_theta, denom_mr real(r_def) :: denom_mr_tmp real(r_def) :: denom_total + integer(i_def), parameter :: n_variables = 5 + real(r_def) :: array_norm(n_variables), array_norm_prev(n_variables) + character(str_def) :: array_names(n_variables) - real(r_def), parameter :: tol = 9.e-1_r_def + real(r_def), parameter :: tol = 3.e-1_r_def + real(r_def), parameter :: indiv_tol = 1.2_r_def integer(i_def) :: n, i, t character(len=4) :: pass_str @@ -212,24 +215,32 @@ module tl_test_timesteps_alg_mod function_space_collection%get_fs(mesh, element_order_h, element_order_v, Wtheta) ) call clone_bundle(mr, diff_mr, nummr) - call invoke( assign_field_random_kernel_type( r_u, 1.0_r_def ) , & - assign_field_random_kernel_type( r_rho, 1.0_r_def ) , & - assign_field_random_kernel_type( r_theta, 1.0_r_def ) , & - assign_field_random_kernel_type( r_exner, 1.0_r_def ) ) + ! Set the 'random data' to be the perturbation + call invoke( name = "copy_fields_to_state", & + setval_X(r_u, u ), & + setval_X(r_exner, exner ), & + setval_X(r_theta, theta), & + setval_X(r_rho, rho ) ) + ! Overwrite the moisture fields with potential temperature do i = 1, nummr - call invoke( assign_field_random_kernel_type( r_mr(i), 1.0_r_def ) ) + call invoke( & + setval_X( r_mr(i), theta ), & + inc_a_times_X( 0.000001, r_mr(i) ), & + setval_X( ls_mr(i), ls_theta ), & + inc_a_times_X( 0.000001, ls_mr(i) )) end do + call moist_dyn_factors_alg( ls_moist_dyn, ls_mr ) - gamma_u = 1.e4_r_def - gamma_theta = 1.e2_r_def - gamma_rho = 1.e0_r_def - gamma_exner = 6.5e-1_r_def + gamma_u = 6.5_r_def + gamma_theta = 6.5_r_def + gamma_rho = 6.5_r_def + gamma_exner = 6.5_r_def - if ( use_moisture ) then - gamma_mr = 1.e-1_r_def + if ( use_moisture ) then + gamma_mr = 6.5_r_def else - gamma_mr = 0.0_r_def + gamma_mr = 0._r_def end if do n = 1, number_gamma_values @@ -277,6 +288,10 @@ module tl_test_timesteps_alg_mod !------------------------------------------------------------------------ ! Nonlinear model run with perturbed state !------------------------------------------------------------------------ + write(log_scratch_space,'(A, I2)') & + "Perturbed Nonlinear, timestep: ", t + call log_event(log_scratch_space, LOG_LEVEL_INFO) + call invoke( setval_X( u, n1_u ), & setval_X( theta, n1_theta ), & setval_X( rho, n1_rho ), & @@ -314,6 +329,10 @@ module tl_test_timesteps_alg_mod !------------------------------------------------------------------------ ! Linear model run !------------------------------------------------------------------------ + write(log_scratch_space,'(A, I2)') & + "Linear, timestep: ", t + call log_event(log_scratch_space, LOG_LEVEL_INFO) + call invoke( setval_X( u, p_u ), & setval_X( theta, p_theta ), & setval_X( rho, p_rho ), & @@ -337,6 +356,7 @@ module tl_test_timesteps_alg_mod call moist_dyn_factors_alg( ls_moist_dyn, ls_mr ) endif + call finalise_linear_model() call initialise_linear_model( mesh, & modeldb ) @@ -359,6 +379,10 @@ module tl_test_timesteps_alg_mod !------------------------------------------------------------------------ ! Nonlinear model run with un-perturbed state !------------------------------------------------------------------------ + write(log_scratch_space,'(A, I2)') & + "Unperturbed NonLinear, timestep: ", t + call log_event(log_scratch_space, LOG_LEVEL_INFO) + call invoke( setval_X( u, n2_u ), & setval_X( theta, n2_theta ), & setval_X( rho, n2_rho ), & @@ -386,6 +410,10 @@ module tl_test_timesteps_alg_mod setval_X( n2_rho, rho ), & setval_X( n2_exner, exner ) ) + do i = 1, nummr + call invoke( setval_X( n2_mr(i), mr(i) ) ) + end do + do i = 1, nummr call invoke( setval_X( ls_mr(i), mr(i) ) ) end do @@ -417,160 +445,67 @@ module tl_test_timesteps_alg_mod norm_mr = norm_mr + norm_mr_tmp end do - norm_diff = norm_rho + norm_theta + norm_exner + norm_u + norm_mr - norm_diff = sqrt( norm_diff ) - norm_u = sqrt( norm_u ) - norm_rho= sqrt( norm_rho ) - norm_theta = sqrt( norm_theta ) - norm_exner = sqrt( norm_exner ) - norm_mr = sqrt( norm_mr ) - - call invoke( X_innerproduct_X( denom_u, p_u ) ) - call invoke( X_innerproduct_X( denom_rho, p_rho ) ) - call invoke( X_innerproduct_X( denom_exner, p_exner ) ) - call invoke( X_innerproduct_X( denom_theta, p_theta ) ) + array_norm(1) = norm_u + array_norm(2) = norm_rho + array_norm(3) = norm_theta + array_norm(4) = norm_mr + array_norm(5) = norm_exner + array_names(1) = 'norm_u' + array_names(2) = 'norm_rho' + array_names(3) = 'norm_theta' + array_names(4) = 'norm_mr' + array_names(5) = 'norm_exner' - denom_mr=0.0_r_def - do i = 1, nummr - call invoke( & - X_innerproduct_X( denom_mr_tmp, p_mr(i) ) ) - denom_mr = denom_mr + denom_mr_tmp - end do + if (n == 2) then + call array_convergence_rate_check( array_norm, array_norm_prev, & + array_names, n_variables, label, tol=tol, indiv_tol=indiv_tol ) + end if - denom_total = denom_rho + denom_theta + denom_exner + denom_u + denom_mr - denom_total = sqrt( denom_total ) - denom_u = sqrt( denom_u ) - denom_rho = sqrt( denom_rho ) - denom_theta = sqrt( denom_theta ) - denom_exner = sqrt( denom_exner ) - denom_mr = sqrt( denom_mr ) + array_norm_prev = array_norm !---------------------------------------------------------------------------- ! Print values out to enable plotting a graph !---------------------------------------------------------------------------- - write( log_scratch_space, & - '( A, E32.12, A, E32.12 )' ) & - ' gamma_total = ' , gamma_u, & - ' norm = ' , norm_diff / denom_total - - call log_event( log_scratch_space, LOG_LEVEL_INFO ) write( log_scratch_space, & '(A, E32.12, A, E32.12)' ) & - ' gamma_rho = ' , gamma_u, & - ' norm = ' , norm_rho / denom_rho + ' gamma_rho = ' , gamma_rho, & + ' norm = ' , norm_rho call log_event( log_scratch_space, LOG_LEVEL_INFO ) write( log_scratch_space, & '(A, E32.12, A, E32.12)' ) & - ' gamma_exner = ' , gamma_u, & - ' norm = ' , norm_exner / denom_exner + ' gamma_exner = ' , gamma_exner, & + ' norm = ' , norm_exner call log_event( log_scratch_space, LOG_LEVEL_INFO ) write( log_scratch_space, & '(A, E32.12, A, E32.12)' ) & - ' gamma_theta = ' , gamma_u, & - ' norm = ' , norm_theta / denom_theta + ' gamma_theta = ' , gamma_theta, & + ' norm = ' , norm_theta call log_event( log_scratch_space, LOG_LEVEL_INFO ) write( log_scratch_space, & '(A, E32.12, A, E32.12)' ) & ' gamma_u = ' , gamma_u, & - ' norm = ' , norm_u / denom_u + ' norm = ' , norm_u call log_event( log_scratch_space, LOG_LEVEL_INFO ) write( log_scratch_space, & '(A, E32.12, A, E32.12)' ) & - ' gamma_mr = ' , gamma_u, & - ' norm = ' , norm_mr / denom_mr + ' gamma_mr = ' , gamma_mr, & + ' norm = ' , norm_mr call log_event( log_scratch_space, LOG_LEVEL_INFO ) -!---------------------------------------------------------------------------- -! Test that gives a PASS or FAIL -!---------------------------------------------------------------------------- - if (n == 2) then - - pass_str = "PASS" - - write( log_scratch_space, '(A)' ) & - "TL Test: " // trim(label) - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - write( log_scratch_space, '(A)' ) "Total " - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - call convergence_pass_string( norm_diff, norm_diff_prev, pass_str, tol ) - - write( log_scratch_space, '(A)' ) "Theta " - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - call convergence_pass_string( norm_theta, norm_theta_prev, pass_str, tol ) - - write( log_scratch_space, '(A)' ) "Rho " - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - call convergence_pass_string( norm_rho, norm_rho_prev, pass_str, tol ) - - write( log_scratch_space, '(A)' ) "Exner " - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - call convergence_pass_string( norm_exner, norm_exner_prev, pass_str, tol ) - - write( log_scratch_space, '(A)' ) "u " - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - call convergence_pass_string( norm_u, norm_u_prev, pass_str, tol ) - - if ( use_moisture ) then - write( log_scratch_space, '(A)' ) "mr " - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - call convergence_pass_string( norm_mr, norm_mr_prev, pass_str, tol ) - end if - - write( log_scratch_space,'(" test",A32," : ",A4)' ) trim(label), pass_str - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - end if - - norm_diff_prev = norm_diff - norm_rho_prev = norm_rho - norm_theta_prev = norm_theta - norm_exner_prev = norm_exner - norm_u_prev = norm_u - norm_mr_prev = norm_mr - end do ! Loop over values of gamma end subroutine test_timesteps - !> @brief Test the convergence rate. - !> @details If the convergence rate is not close to 4, within the - !! specified tolerance, then change the pass string to FAIL. - !> @param[in] norm_diff Current norm - !> @param[in] norm_diff_prev Previous norm - !> @param[inout] pass_str Pass string (either PASS or FAIL) - !> @param[in] tol Test Tolerance - subroutine convergence_pass_string( norm_diff, norm_diff_prev, pass_str, tol ) - implicit none - - real(r_def), intent(in) :: norm_diff, norm_diff_prev - character(len=4), intent(inout) :: pass_str - real(r_def), intent(in) :: tol - - real(r_def) :: conv_rate - conv_rate = norm_diff_prev / norm_diff - - conv_rate = norm_diff_prev / norm_diff - write( log_scratch_space, '(A, E16.8)') & - "Convergence rate: ", conv_rate - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - if ( abs( conv_rate - 4.0_r_def ) >= tol ) then - pass_str = "FAIL" - endif - - end subroutine convergence_pass_string - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Converts a string to a timestep number. !> diff --git a/science/linear/integration-test/tl_test/tl_test_timesteps_random_alg_mod.x90 b/science/linear/integration-test/tl_test/tl_test_timesteps_random_alg_mod.x90 new file mode 100644 index 00000000..9b8ed013 --- /dev/null +++ b/science/linear/integration-test/tl_test/tl_test_timesteps_random_alg_mod.x90 @@ -0,0 +1,595 @@ +!----------------------------------------------------------------------------- +! (C) Crown copyright 2026 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!----------------------------------------------------------------------------- + +!>@brief The tangent linear test for multiple timesteps +!! (Taylor remainder convergence). +!>@details Test whether the tangent linear code is tangent linear +!! to the corresponding nonlinear code, for multiple timesteps, +!! using the Taylor remainder convergence test. +module tl_test_timesteps_random_alg_mod + use sci_assign_field_random_kernel_mod, & + only: assign_field_random_kernel_type + use constants_mod, only: i_timestep, i_def, r_def, l_def, & + str_def + use field_array_mod, only: field_array_type + use field_mod, only: field_type + use sci_field_bundle_builtins_mod, only: clone_bundle + use field_collection_mod, only: field_collection_type + use finite_element_config_mod, only: element_order_h, element_order_v + use formulation_config_mod, only: moisture_formulation, & + moisture_formulation_dry + use fs_continuity_mod, only: W2, W3, Wtheta + use function_space_collection_mod, only: function_space_collection + use gungho_step_mod, only: gungho_step + use gungho_model_mod, only: initialise_model, & + finalise_model + use driver_modeldb_mod, only: modeldb_type + use log_mod, only: log_event, & + log_scratch_space, & + LOG_LEVEL_INFO, & + LOG_LEVEL_ERROR + use linear_step_mod, only: linear_step + use linear_model_mod, only: initialise_linear_model, & + finalise_linear_model + use mesh_mod, only: mesh_type + use model_clock_mod, only: model_clock_type + use moist_dyn_mod, only: num_moist_factors + use mr_indices_mod, only: nummr + use semi_implicit_solver_alg_mod, only: semi_implicit_solver_alg_final + use si_operators_alg_mod, only: final_si_operators + use timestepping_config_mod, only: dt, & + method, & + method_rk + use time_config_mod, only: timestep_start, timestep_end + use tl_test_convergence_rate_check, only: convergence_rate_check + use validity_test_config_mod, only: number_gamma_values, & + update_ls_frequency + use moist_dyn_factors_alg_mod, only: moist_dyn_factors_alg + use tl_moist_dyn_factors_alg_mod, only: tl_moist_dyn_factors_alg + + implicit none + + private convergence_pass_string + public test_timesteps_random + + contains + + !> @brief Test the tangent linear runge-kutta algorithm. + !> @param[in] modeldb The working data set for a model run + !> @param[in] mesh The current 3d mesh + !> @param[in] twod_mesh The current 2d mesh + !> @param[in] model_clock Time within the model. + !> + subroutine test_timesteps_random( modeldb, & + mesh, & + twod_mesh, & + model_clock ) + + implicit none + + type(modeldb_type), target, intent(inout) :: modeldb + type(mesh_type), pointer, intent(in) :: mesh + type(mesh_type), pointer, intent(in) :: twod_mesh + class(model_clock_type), intent(in) :: model_clock + + character(str_def) :: label = "timesteps" + + type(field_collection_type), pointer :: prognostic_fields + type(field_collection_type), pointer :: ls_fields + + type(field_type), pointer :: ls_u + type(field_type), pointer :: ls_rho + type(field_type), pointer :: ls_theta + type(field_type), pointer :: ls_exner + type(field_type), pointer :: ls_moist_dyn(:) + type(field_type), pointer :: ls_mr(:) + + type(field_type), pointer :: u + type(field_type), pointer :: rho + type(field_type), pointer :: theta + type(field_type), pointer :: exner + type(field_type), pointer :: moist_dyn(:) + type(field_type), pointer :: mr(:) + + type(field_type), dimension(nummr) :: p_mr + type(field_type), dimension(nummr) :: n1_mr + type(field_type), dimension(nummr) :: n2_mr + type(field_type), dimension(nummr) :: r_mr + type(field_type), dimension(nummr) :: diff_mr + type(field_type) :: p_u, p_rho, p_theta, p_exner + type(field_type) :: n1_u, n1_rho, n1_theta, n1_exner + type(field_type) :: n2_u, n2_rho, n2_theta, n2_exner + type(field_type) :: r_u, r_rho, r_theta, r_exner + type(field_type) :: diff_u, diff_rho, diff_theta, diff_exner + + type(field_collection_type), pointer :: moisture_fields + type(field_array_type), pointer :: mr_array + type(field_array_type), pointer :: moist_dyn_array + type(field_array_type), pointer :: ls_mr_array + type(field_array_type), pointer :: ls_moist_dyn_array + + real(r_def) :: gamma_u, gamma_rho, gamma_exner, gamma_theta, gamma_mr + real(r_def) :: norm_u, norm_rho, norm_exner, norm_theta, norm_mr + real(r_def) :: norm_mr_tmp + real(r_def) :: norm_diff, norm_diff_prev + real(r_def) :: norm_u_prev, norm_rho_prev, norm_exner_prev + real(r_def) :: norm_theta_prev, norm_mr_prev + + real(r_def), parameter :: tol = 9.e-1_r_def + + integer(i_def) :: n, i, t + character(len=4) :: pass_str + logical(l_def) :: use_moisture + + call log_event( "TL Test: " // trim(label), & + LOG_LEVEL_INFO ) + + use_moisture = ( moisture_formulation /= moisture_formulation_dry ) + + if ( method == method_rk .and. use_moisture ) then + call log_event( & + 'Not possible to use Runge Kutta time stepping with moisture', & + log_level_error ) + end if + + prognostic_fields => modeldb%fields%get_field_collection( & + "prognostic_fields") + ls_fields => modeldb%fields%get_field_collection("ls_fields") + moisture_fields => modeldb%fields%get_field_collection("moisture_fields") + call moisture_fields%get_field("mr", mr_array) + call moisture_fields%get_field("moist_dyn", moist_dyn_array) + mr => mr_array%bundle + moist_dyn => moist_dyn_array%bundle + call moisture_fields%get_field("ls_mr", ls_mr_array) + call moisture_fields%get_field("ls_moist_dyn", ls_moist_dyn_array) + ls_mr => ls_mr_array%bundle + ls_moist_dyn => ls_moist_dyn_array%bundle + + ! Input + call ls_fields%get_field('ls_u', ls_u) + call ls_fields%get_field('ls_rho', ls_rho) + call ls_fields%get_field('ls_theta', ls_theta) + call ls_fields%get_field('ls_exner', ls_exner) + + call prognostic_fields%get_field('u', u) + call prognostic_fields%get_field('rho', rho) + call prognostic_fields%get_field('theta', theta) + call prognostic_fields%get_field('exner', exner) + + call r_u%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W2) ) + call r_rho%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call r_exner%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call r_theta%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, Wtheta) ) + call clone_bundle(mr, r_mr, nummr) + + call p_u%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W2) ) + call p_rho%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call p_exner%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call p_theta%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, Wtheta) ) + call clone_bundle(mr, p_mr, nummr) + + call n1_u%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W2) ) + call n1_rho%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call n1_exner%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call n1_theta%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, Wtheta) ) + call clone_bundle(mr, n1_mr, nummr) + + call n2_u%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W2) ) + call n2_rho%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call n2_exner%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call n2_theta%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, Wtheta) ) + call clone_bundle(mr, n2_mr, nummr) + + call diff_u%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W2) ) + call diff_rho%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call diff_exner%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, W3) ) + call diff_theta%initialise( vector_space = & + function_space_collection%get_fs(mesh, element_order_h, element_order_v, Wtheta) ) + call clone_bundle(mr, diff_mr, nummr) + + call invoke( assign_field_random_kernel_type( r_u, 1.0_r_def ) , & + assign_field_random_kernel_type( r_rho, 1.0_r_def ) , & + assign_field_random_kernel_type( r_theta, 1.0_r_def ) , & + assign_field_random_kernel_type( r_exner, 1.0_r_def ) ) + + do i = 1, nummr + call invoke( assign_field_random_kernel_type( r_mr(i), 1.0_r_def ) ) + end do + + gamma_u = 1.e4_r_def + gamma_theta = 1.0_r_def + gamma_rho = 1.e0_r_def + gamma_exner = 6.5e-1_r_def + + if ( use_moisture ) then + gamma_mr = 1.e-2_r_def + else + gamma_mr = 0.0_r_def + end if + + gamma_u = gamma_u /24.0_r_def + gamma_theta = gamma_theta / 24.0_r_def + gamma_rho = gamma_rho / 24.0_r_def + gamma_exner = gamma_exner / 24.0_r_def + gamma_mr = gamma_mr / 24.0_r_def + + do n = 1, number_gamma_values + gamma_u = gamma_u / 2.0_r_def + gamma_theta = gamma_theta / 2.0_r_def + gamma_rho = gamma_rho / 2.0_r_def + gamma_exner = gamma_exner / 2.0_r_def + gamma_mr = gamma_mr / 2.0_r_def + + ! Initial conditions for nonlinear perturbed + call invoke( aX_plus_Y( n1_u, gamma_u, r_u, ls_u ), & + aX_plus_Y( n1_theta, gamma_theta, r_theta, ls_theta ), & + aX_plus_Y( n1_rho, gamma_rho, r_rho, ls_rho ), & + aX_plus_Y( n1_exner, gamma_exner, r_exner, ls_exner ) ) + + do i = 1, nummr + call invoke( aX_plus_Y( n1_mr(i), gamma_mr, & + r_mr(i), ls_mr(i) ) ) + end do + + ! Initial conditions for nonlinear unperturbed + call invoke( setval_X( n2_u, ls_u ), & + setval_X( n2_theta, ls_theta ), & + setval_X( n2_rho, ls_rho ), & + setval_X( n2_exner, ls_exner ) ) + + do i = 1, nummr + call invoke( setval_X( n2_mr(i), ls_mr(i) ) ) + end do + + ! Initial conditions for the linear + call invoke( a_times_X( p_u, gamma_u, r_u ), & + a_times_X( p_theta, gamma_theta, r_theta ), & + a_times_X( p_rho, gamma_rho, r_rho ), & + a_times_X( p_exner, gamma_exner, r_exner ) ) + + do i = 1, nummr + call invoke( a_times_X( p_mr(i), gamma_mr, & + r_mr(i) ) ) + end do + + ! Multiple timesteps + do t = parse_instance(timestep_start), parse_instance(timestep_end) + +!------------------------------------------------------------------------ +! Nonlinear model run with perturbed state +!------------------------------------------------------------------------ + call invoke( setval_X( u, n1_u ), & + setval_X( theta, n1_theta ), & + setval_X( rho, n1_rho ), & + setval_X( exner, n1_exner ) ) + do i = 1, nummr + call invoke( setval_X( mr(i), n1_mr(i) ) ) + enddo + call moist_dyn_factors_alg( moist_dyn, mr ) + + ! Clean up solvers and timestep method left over from running + ! the model to set up model state and linear state + call final_si_operators() + call semi_implicit_solver_alg_final + call finalise_model(modeldb) + + call initialise_model( mesh, & + modeldb ) + + call gungho_step( mesh, & + twod_mesh, & + modeldb, & + model_clock ) + + call finalise_model(modeldb) + + call invoke( setval_X( n1_u, u ), & + setval_X( n1_theta, theta ), & + setval_X( n1_rho, rho ), & + setval_X( n1_exner, exner) ) + + do i = 1, nummr + call invoke( setval_X( n1_mr(i), mr(i) ) ) + enddo + +!------------------------------------------------------------------------ +! Linear model run +!------------------------------------------------------------------------ + call invoke( setval_X( u, p_u ), & + setval_X( theta, p_theta ), & + setval_X( rho, p_rho ), & + setval_X( exner, p_exner) ) + do i = 1, nummr + call invoke( setval_X( mr(i), p_mr(i) ) ) + enddo + call tl_moist_dyn_factors_alg( moist_dyn, mr ) + + ! Set the linearisation state + ! The linearisation state is copied from the nonlinear run, + ! every few timesteps. + if ( mod(t, update_ls_frequency) == 0 ) then + call invoke( setval_X( ls_u, n2_u ), & + setval_X( ls_theta, n2_theta ), & + setval_X( ls_rho, n2_rho ), & + setval_X( ls_exner, n2_exner) ) + do i = 1, nummr + call invoke( setval_X( ls_mr(i), n2_mr(i) ) ) + enddo + call moist_dyn_factors_alg( ls_moist_dyn, ls_mr ) + endif + + call finalise_linear_model() + call initialise_linear_model( mesh, & + modeldb ) + + call linear_step( mesh, & + twod_mesh, & + modeldb, & + model_clock ) + + call finalise_linear_model() + + + call invoke( setval_X( p_u, u ), & + setval_X( p_theta, theta ), & + setval_X( p_rho, rho ), & + setval_X( p_exner, exner) ) + do i = 1, nummr + call invoke( setval_X( p_mr(i), mr(i) ) ) + enddo + +!------------------------------------------------------------------------ +! Nonlinear model run with un-perturbed state +!------------------------------------------------------------------------ + call invoke( setval_X( u, n2_u ), & + setval_X( theta, n2_theta ), & + setval_X( rho, n2_rho ), & + setval_X( exner, n2_exner ) ) + + do i = 1, nummr + call invoke( setval_X( mr(i), n2_mr(i) ) ) + end do + + call moist_dyn_factors_alg( moist_dyn, mr ) + + + call initialise_model( mesh, & + modeldb ) + + call gungho_step( mesh, & + twod_mesh, & + modeldb, & + model_clock ) + + call finalise_model(modeldb) + + call invoke( setval_X( n2_u, u ), & + setval_X( n2_theta, theta ), & + setval_X( n2_rho, rho ), & + setval_X( n2_exner, exner ) ) + + do i = 1, nummr + call invoke( setval_X( n2_mr(i), mr(i) ) ) + call invoke( setval_X( ls_mr(i), mr(i) ) ) + end do + + end do ! timesteps loop + +!------------------------------------------------------------------------ +! Calculate norms +!------------------------------------------------------------------------ + call invoke( X_minus_Y( diff_u, n1_u, n2_u ), & + inc_X_minus_Y( diff_u, p_u ), & + X_innerproduct_X( norm_u, diff_u ) ) + call invoke( X_minus_Y( diff_rho, n1_rho, n2_rho ), & + inc_X_minus_Y( diff_rho, p_rho ), & + X_innerproduct_X( norm_rho, diff_rho ) ) + call invoke( X_minus_Y( diff_exner, n1_exner, n2_exner ), & + inc_X_minus_Y( diff_exner, p_exner ), & + X_innerproduct_X( norm_exner, diff_exner ) ) + call invoke( X_minus_Y( diff_theta, n1_theta, n2_theta ), & + inc_X_minus_Y( diff_theta, p_theta ), & + X_innerproduct_X( norm_theta, diff_theta ) ) + + norm_mr=0.0_r_def + do i = 1, nummr + call invoke( & + X_minus_Y( diff_mr(i), n1_mr(i), n2_mr(i) ), & + inc_X_minus_Y( diff_mr(i), p_mr(i) ), & + X_innerproduct_X( norm_mr_tmp, diff_mr(i) ) ) + norm_mr = norm_mr + norm_mr_tmp + end do + + norm_diff = norm_rho + norm_theta + norm_exner + norm_u + norm_mr + +!---------------------------------------------------------------------------- +! Print values out to enable plotting a graph +!---------------------------------------------------------------------------- + write( log_scratch_space, & + '( A, E32.12, A, E32.12 )' ) & + ' gamma_total = ' , gamma_u, & + ' norm = ' , norm_diff + + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + write( log_scratch_space, & + '(A, E32.12, A, E32.12)' ) & + ' gamma_rho = ' , gamma_u, & + ' norm = ' , norm_rho + + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + write( log_scratch_space, & + '(A, E32.12, A, E32.12)' ) & + ' gamma_exner = ' , gamma_u, & + ' norm = ' , norm_exner + + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + write( log_scratch_space, & + '(A, E32.12, A, E32.12)' ) & + ' gamma_theta = ' , gamma_u, & + ' norm = ' , norm_theta + + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + write( log_scratch_space, & + '(A, E32.12, A, E32.12)' ) & + ' gamma_u = ' , gamma_u, & + ' norm = ' , norm_u + + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + write( log_scratch_space, & + '(A, E32.12, A, E32.12)' ) & + ' gamma_mr = ' , gamma_u, & + ' norm = ' , norm_mr + + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + +!---------------------------------------------------------------------------- +! Test that gives a PASS or FAIL +!---------------------------------------------------------------------------- + if (n == 2) then + + pass_str = "PASS" + + write( log_scratch_space, '(A)' ) & + "TL Test: " // trim(label) + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + write( log_scratch_space, '(A)' ) "Total " + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + call convergence_pass_string( norm_diff, norm_diff_prev, pass_str, tol ) + + write( log_scratch_space, '(A)' ) "Theta " + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + call convergence_pass_string( norm_theta, norm_theta_prev, pass_str, tol ) + + write( log_scratch_space, '(A)' ) "Rho " + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + call convergence_pass_string( norm_rho, norm_rho_prev, pass_str, tol ) + + write( log_scratch_space, '(A)' ) "Exner " + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + call convergence_pass_string( norm_exner, norm_exner_prev, pass_str, tol ) + + write( log_scratch_space, '(A)' ) "u " + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + call convergence_pass_string( norm_u, norm_u_prev, pass_str, tol ) + + if ( use_moisture ) then + write( log_scratch_space, '(A)' ) "mr " + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + call convergence_pass_string( norm_mr, norm_mr_prev, pass_str, tol ) + end if + + write( log_scratch_space,'(" test",A32," : ",A4)' ) trim(label), pass_str + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + end if + + norm_diff_prev = norm_diff + norm_rho_prev = norm_rho + norm_theta_prev = norm_theta + norm_exner_prev = norm_exner + norm_u_prev = norm_u + norm_mr_prev = norm_mr + + end do ! Loop over values of gamma + + end subroutine test_timesteps_random + + !> @brief Test the convergence rate. + !> @details If the convergence rate is not close to 4, within the + !! specified tolerance, then change the pass string to FAIL. + !> @param[in] norm_diff Current norm + !> @param[in] norm_diff_prev Previous norm + !> @param[inout] pass_str Pass string (either PASS or FAIL) + !> @param[in] tol Test Tolerance + subroutine convergence_pass_string( norm_diff, norm_diff_prev, pass_str, tol ) + implicit none + + real(r_def), intent(in) :: norm_diff, norm_diff_prev + character(len=4), intent(inout) :: pass_str + real(r_def), intent(in) :: tol + + real(r_def) :: conv_rate + conv_rate = norm_diff_prev / norm_diff + + conv_rate = norm_diff_prev / norm_diff + write( log_scratch_space, '(A, E16.8)') & + "Convergence rate: ", conv_rate + call log_event( log_scratch_space, LOG_LEVEL_INFO ) + + if ( abs( conv_rate - 4.0_r_def ) >= tol ) then + pass_str = "FAIL" + endif + + end subroutine convergence_pass_string + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Converts a string to a timestep number. + !> + !> @param[in] string Correctly formatted timestep number. + !> @return Timestep number. + !> + function parse_instance( string ) result(instance) + + implicit none + + character(*), intent(in) :: string + integer(i_timestep) :: instance + + integer :: string_size + integer :: format_size + integer :: status + character(:), allocatable :: fmt + + ! The "I0" formatting string does not work for input. Instead the exact + ! number of digits to be read must be requested. Thus we need to construct + ! a format string from the size of the string. + ! + string_size = len(string) + format_size = string_size + 3 + allocate( character(format_size) :: fmt, stat=status ) + if ( status /= 0 ) then + call log_event( & + 'TL test parse_instance: Unable to allocate format', & + log_level_error ) + end if + write( fmt, '("(I", I0, ")")' ) string_size + read( string, fmt ) instance + deallocate( fmt ) + + if ( instance < 0 ) then + call log_event( & + 'TL test parse_instance: Instances may not be negative', & + log_level_error ) + end if + + end function parse_instance + +end module tl_test_timesteps_random_alg_mod diff --git a/science/linear/integration-test/tl_test/tl_test_transport_control_mod.x90 b/science/linear/integration-test/tl_test/tl_test_transport_control_mod.x90 index 02d96d2e..d74fd3df 100644 --- a/science/linear/integration-test/tl_test/tl_test_transport_control_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_transport_control_mod.x90 @@ -21,7 +21,6 @@ module tl_test_transport_control_mod use function_space_collection_mod, only: function_space_collection use derived_config_mod, only: bundle_size use driver_modeldb_mod, only: modeldb_type - use moist_dyn_mod, only: num_moist_factors use mr_indices_mod, only: nummr use log_mod, only: log_event, & log_scratch_space, & @@ -31,7 +30,7 @@ module tl_test_transport_control_mod only: gungho_transport_control_alg_init, & gungho_transport_control_alg use tl_transport_control_alg_mod, only: tl_transport_control_alg - use tl_test_convergence_rate_check, only: convergence_rate_check + use tl_test_convergence_rate_check, only: array_convergence_rate_check implicit none @@ -61,10 +60,11 @@ module tl_test_transport_control_mod type(field_type), pointer :: ls_u => null() type(field_type), pointer :: ls_rho => null() type(field_type), pointer :: ls_theta => null() - type(field_type), pointer :: ls_exner => null() - type(field_type), pointer :: ls_moist_dyn(:) => null() type(field_type), pointer :: ls_mr(:) => null() - type(field_type), pointer :: moist_dyn(:) => null() + type(field_type), pointer :: ls_exner => null() + type(field_type), pointer :: u => null() + type(field_type), pointer :: rho => null() + type(field_type), pointer :: theta => null() type(field_type), pointer :: mr(:) => null() type(field_type) :: state(bundle_size) @@ -75,10 +75,7 @@ module tl_test_transport_control_mod type(field_type) :: n1_rhs(bundle_size) type(field_type) :: n2_rhs(bundle_size) type(field_type) :: diff(bundle_size) - type(field_type), dimension(num_moist_factors) :: p_moist_dyn - type(field_type), dimension(num_moist_factors) :: n_moist_dyn - type(field_type), dimension(num_moist_factors) :: r_moist_dyn - type(field_type), dimension(num_moist_factors) :: diff_moist_dyn + type(field_type), dimension(nummr) :: p_mr_in type(field_type), dimension(nummr) :: n_mr_in type(field_type), dimension(nummr) :: ls_mr_out @@ -91,36 +88,32 @@ module tl_test_transport_control_mod type(field_collection_type), pointer :: moisture_fields => null() type(field_array_type), pointer :: mr_array => null() - type(field_array_type), pointer :: moist_dyn_array => null() type(field_array_type), pointer :: ls_mr_array => null() - type(field_array_type), pointer :: ls_moist_dyn_array => null() - - real(r_def) :: gamma_u, gamma_rho, gamma_exner, gamma_theta - real(r_def) :: gamma_moist_dyn, gamma_mr - real(r_def) :: norm_u, norm_rho, norm_exner, norm_theta - real(r_def) :: norm_moist_dyn, norm_mr - real(r_def) :: norm_diff, norm_diff_prev - real(r_def) :: norm_moist_dyn_tmp, norm_mr_tmp - real(r_def), parameter :: tol = 1.e-2_r_def + real(r_def) :: gamma_u, gamma_rho, gamma_theta + real(r_def) :: gamma_mr + real(r_def) :: norm_u, norm_rho, norm_theta + real(r_def) :: norm_mr + real(r_def) :: norm_mr_tmp + integer(i_def), parameter :: n_variables = 4 + real(r_def) :: array_norm(n_variables), array_norm_prev(n_variables) + character(str_def) :: array_names(n_variables) + real(r_def), parameter :: tol = 5.e-2_r_def + real(r_def), parameter :: indiv_tol = 5.e-1_r_def integer :: n, outer, i call log_event( "TL Test: " // trim(label), & - LOG_LEVEL_INFO ) + LOG_LEVEL_INFO ) prognostic_fields => modeldb%fields%get_field_collection( & "prognostic_fields") ls_fields => modeldb%fields%get_field_collection("ls_fields") moisture_fields => modeldb%fields%get_field_collection("moisture_fields") call moisture_fields%get_field("mr", mr_array) - call moisture_fields%get_field("moist_dyn", moist_dyn_array) mr => mr_array%bundle - moist_dyn => moist_dyn_array%bundle call moisture_fields%get_field("ls_mr", ls_mr_array) - call moisture_fields%get_field("ls_moist_dyn", ls_moist_dyn_array) ls_mr => ls_mr_array%bundle - ls_moist_dyn => ls_moist_dyn_array%bundle ! Input call ls_fields%get_field('ls_u', ls_u) @@ -128,6 +121,12 @@ module tl_test_transport_control_mod call ls_fields%get_field('ls_theta', ls_theta) call ls_fields%get_field('ls_exner', ls_exner) + ! Overwrite the moisture fields with potential temperature + do i = 1, nummr + call invoke( & + setval_X( ls_mr(i), ls_theta)) + end do + call ls_state(igh_u)%initialise( vector_space = ls_u%get_function_space() ) call ls_state(igh_t)%initialise( vector_space = ls_theta%get_function_space() ) call ls_state(igh_d)%initialise( vector_space = ls_rho%get_function_space() ) @@ -144,11 +143,6 @@ module tl_test_transport_control_mod call clone_bundle(ls_state, n2_rhs, bundle_size) call clone_bundle(ls_state, diff, bundle_size) - call clone_bundle(moist_dyn, p_moist_dyn, num_moist_factors) - call clone_bundle(moist_dyn, r_moist_dyn, num_moist_factors) - call clone_bundle(moist_dyn, n_moist_dyn, num_moist_factors) - call clone_bundle(moist_dyn, diff_moist_dyn, num_moist_factors) - call clone_bundle(mr, r_mr, nummr) call clone_bundle(mr, p_mr_in, nummr) call clone_bundle(mr, n_mr_in, nummr) @@ -161,23 +155,36 @@ module tl_test_transport_control_mod setval_X(ls_state(igh_u), ls_u ), & setval_X(ls_state(igh_t), ls_theta), & setval_X(ls_state(igh_d), ls_rho ), & - setval_X(ls_state(igh_p), ls_exner), & + setval_X(ls_state(igh_p), ls_exner ), & setval_X(ls_advected_u, ls_u ) ) call gungho_transport_control_alg_init( mesh ) - call invoke( assign_field_random_kernel_type( random(igh_u), 1.0_r_def ) , & - assign_field_random_kernel_type( random(igh_d), 1.0_r_def ) , & - assign_field_random_kernel_type( random(igh_t), 1.0_r_def ) , & - assign_field_random_kernel_type( random(igh_p), 1.0_r_def ) ) + ! Perturbation + call prognostic_fields%get_field('u', u) + call prognostic_fields%get_field('rho', rho) + call prognostic_fields%get_field('theta', theta) + + ! Set the 'random perturbation as the ls_state plus pert' + call invoke( name = "copy_fields_to_state", & + setval_X(random(igh_u), ls_u ), & + setval_X(random(igh_t), ls_theta), & + setval_X(random(igh_d), ls_rho ) ) + + do i = 1, nummr + call invoke( & + setval_X( r_mr(i), ls_mr(i) )) + end do - do i = 1, num_moist_factors - call invoke( assign_field_random_kernel_type( r_moist_dyn(i), 1.0_r_def ) ) - enddo + call invoke( name = "inc", & + inc_X_plus_Y(random(igh_u), u ), & + inc_X_plus_Y(random(igh_t), theta), & + inc_X_plus_Y(random(igh_d), rho ) ) do i = 1, nummr - call invoke( assign_field_random_kernel_type( r_mr(i), 1.0_r_def ) ) - enddo + call invoke( & + inc_X_plus_Y( r_mr(i), theta )) + end do call set_bundle_scalar( 0.0_r_def, n1_rhs, bundle_size ) @@ -192,19 +199,15 @@ module tl_test_transport_control_mod outer, & cheap_update = .false.) - gamma_u = 7.e8_r_def - gamma_theta = 1.e5_r_def - gamma_rho = 1.e5_r_def - gamma_exner = 1.e5_r_def - gamma_moist_dyn=1.e0_r_def - gamma_mr = 1.e0_r_def + gamma_u = 1._r_def + gamma_theta = 1._r_def + gamma_rho = 1._r_def + gamma_mr = 1._r_def do n=1,2 gamma_u = gamma_u/2.0_r_def gamma_theta = gamma_theta/2.0_r_def gamma_rho = gamma_rho/2.0_r_def - gamma_exner = gamma_exner/2.0_r_def - gamma_moist_dyn = gamma_moist_dyn/2.0_r_def gamma_mr = gamma_mr/2.0_r_def call set_bundle_scalar( 0.0_r_def, n2_rhs, bundle_size ) @@ -212,7 +215,7 @@ module tl_test_transport_control_mod call invoke( a_times_X( p_state(igh_u), gamma_u, random(igh_u) ), & a_times_X( p_state(igh_t), gamma_theta, random(igh_t) ), & a_times_X( p_state(igh_d), gamma_rho, random(igh_d) ), & - a_times_X( p_state(igh_p), gamma_exner, random(igh_p) ), & + setval_c( p_state(igh_p), 0.0_r_def ), & a_times_X( p_advected_u, gamma_u, random(igh_u) ) ) ! state = ls_state + p_state @@ -221,13 +224,6 @@ module tl_test_transport_control_mod call invoke( setval_X( n_advected_u, ls_advected_u), & inc_X_plus_Y(n_advected_u, p_advected_u ) ) - do i = 1, num_moist_factors - call invoke( & - a_times_X( p_moist_dyn(i), gamma_moist_dyn, r_moist_dyn(i) ), & - setval_X( n_moist_dyn(i), ls_moist_dyn(i) ), & - inc_X_plus_Y( n_moist_dyn(i), p_moist_dyn(i) ) ) - end do - do i = 1, nummr call invoke( & a_times_X( p_mr_in(i), gamma_mr, r_mr(i) ), & @@ -263,63 +259,43 @@ module tl_test_transport_control_mod call invoke( & inc_X_minus_Y( diff(igh_u), p_rhs(igh_u) ), & inc_X_minus_Y( diff(igh_t), p_rhs(igh_t) ), & - inc_X_minus_Y( diff(igh_d), p_rhs(igh_d) ), & - inc_X_minus_Y( diff(igh_p), p_rhs(igh_p) ) ) + inc_X_minus_Y( diff(igh_d), p_rhs(igh_d) ) ) - call invoke( X_innerproduct_X( norm_u, diff(igh_u) ) , & + call invoke( & + X_innerproduct_X( norm_u, diff(igh_u) ) , & X_innerproduct_X( norm_rho, diff(igh_d) ) , & - X_innerproduct_X( norm_theta, diff(igh_t) ) , & - X_innerproduct_X( norm_exner, diff(igh_p) ) ) + X_innerproduct_X( norm_theta, diff(igh_t) ) ) - norm_moist_dyn=0.0_r_def - do i = 1, num_moist_factors + + do i = 1, nummr call invoke( & - X_minus_Y( diff_moist_dyn(i), n_moist_dyn(i), moist_dyn(i) ), & - inc_X_minus_Y( diff_moist_dyn(i), p_moist_dyn(i) ), & - X_innerproduct_X( norm_moist_dyn_tmp, diff_moist_dyn(i) ) ) - norm_moist_dyn = norm_moist_dyn + norm_moist_dyn_tmp + X_minus_Y( diff_mr(i), n_mr_out(i), ls_mr_out(i) ), & + inc_X_minus_Y( diff_mr(i), p_mr_out(i) ) ) end do norm_mr=0.0_r_def do i = 1, nummr - call invoke( & - X_minus_Y( diff_mr(i), n_mr_out(i), ls_mr_out(i) ), & - inc_X_minus_Y( diff_mr(i), p_mr_out(i) ), & + call invoke( & X_innerproduct_X( norm_mr_tmp, diff_mr(i) ) ) norm_mr = norm_mr + norm_mr_tmp end do - write( log_scratch_space, & - '(A, E32.12, A, E32.12 )' ) & - ' norm_u = ' , norm_u, & - ' norm_rho = ' , norm_rho - - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - write( log_scratch_space, & - '(A, E32.12, A, E32.12 )' ) & - ' norm_theta = ' , norm_theta, & - ' norm_moist_dyn = ' , norm_moist_dyn - - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - write( log_scratch_space, & - '(A, E32.12, A, E32.12 )' ) & - ' norm_exner = ' , norm_exner, & - ' norm_mr = ' , norm_mr - - call log_event( log_scratch_space, LOG_LEVEL_INFO ) - - norm_diff = norm_rho + norm_theta + norm_exner + norm_u & - + norm_moist_dyn + norm_mr - norm_diff = sqrt(norm_diff) + array_norm(1) = norm_u + array_norm(2) = norm_rho + array_norm(3) = norm_theta + array_norm(4) = norm_mr + array_names(1) = 'norm_u' + array_names(2) = 'norm_rho' + array_names(3) = 'norm_theta' + array_names(4) = 'norm_mr' if (n == 2) then - call convergence_rate_check( norm_diff, norm_diff_prev, & - label, tol=tol ) + call array_convergence_rate_check( array_norm, array_norm_prev, & + array_names, n_variables, label, tol=tol, indiv_tol=indiv_tol) end if - norm_diff_prev = norm_diff + array_norm_prev = array_norm + end do end subroutine test_transport_control diff --git a/science/linear/integration-test/tl_test/tl_test_vorticity_mod.x90 b/science/linear/integration-test/tl_test/tl_test_vorticity_mod.x90 index 47e5cdfc..c129dd3f 100644 --- a/science/linear/integration-test/tl_test/tl_test_vorticity_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_vorticity_mod.x90 @@ -203,8 +203,6 @@ module tl_test_vorticity_mod inc_X_minus_Y( diff, p_rhs_u ), & X_innerproduct_X( norm_diff, diff ) ) - norm_diff = sqrt(norm_diff) - write( log_scratch_space, '(A, E32.12, A, E32.12)' ) & 'gamma = ', gamma, ' norm = ', norm_diff diff --git a/science/linear/plot_convergence/plot_convergence.py b/science/linear/plot_convergence/plot_convergence.py new file mode 100644 index 00000000..e5c7412d --- /dev/null +++ b/science/linear/plot_convergence/plot_convergence.py @@ -0,0 +1,120 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +############################################################################## +# (c) Crown copyright 2021 Met Office. All rights reserved. +# The file LICENCE, distributed with this code, contains details of the terms +# under which the code may be used. +############################################################################## +''' +Using the output derived from the test_timesteps integration test with +10 values of gamma, plot a graph of the relative error (linearisation error) +against the size of the perturbation for different prognostic variables +(gamma). +''' +import os +import pandas as pd +import matplotlib.pyplot as plt +import numpy as np + + +def plot_data(filename, axes, variable, color, shape): + ''' + Create a data frame with columns gamma (size of the perturbation) + and norm (the relative error). Plot this data with a log-log axis. + ''' + + datafile = open(filename, 'r') + + norm_line = [] + line = datafile.readline + + while line: + line = datafile.readline() + + if variable in line: + split_line = line.replace('norm', '').replace('\n', '').split('=') + norm_line.append([float(split_line[1]), float(split_line[2])]) + + norm_df = pd.DataFrame(norm_line, columns=['gamma', 'norm']) + + datafile.close() + + # Square root (as the read data is only the inner product and + # has not included the square root to give the norm) + norm_df['norm'] = np.sqrt(norm_df['norm']) + + # Normalise - to give a relative error + normalise = norm_df['norm'].iloc[4] + norm_df['norm'] = norm_df['norm'] / normalise + + # Plot the data + if (CONFIG == 'nwp_gal9'): + ymin = 10**-2 + ymax = 10**2 + xmin = 10**-3 + xmax = 10**1 + elif (CONFIG == 'semi_implicit'): + ymin = 10**-2 + ymax = 10**2 + xmin = 10**-1 + xmax = 10**4 + elif (CONFIG == 'runge_kutta'): + ymin = 10**-3 + ymax = 10**2 + xmin = 10**-1 + xmax = 10**4 + else: + print(CONFIG+' not listed') + + norm_df.plot.scatter(x='gamma', y='norm', loglog=True, + xlim=(xmin, xmax), + ylim=(ymin, ymax), + ax=axes, color=color, + marker=shape, label = variable) + + # Check extremes + norm_min = norm_df['norm'].min() + norm_max = norm_df['norm'].max() + if (norm_min < ymin) : + print('Warning: Min value is'+ str(norm_min)) + if (norm_max > ymax) : + print('Warning: Max value is'+ str(norm_max)) + + # Plot the expected line + centre = norm_df['gamma'].iloc[4] + expected_x = [10**-2 *centre, centre, 100* centre] + expected_y = [10**-2, 10**0, 100] + plt.plot(expected_x, expected_y) + + +def make_plot(directory, filename): + ''' + Plot the data for the different prognostic variables, together with the + expected gradient ( which is linear ), on the same plot. + ''' + + axs = plt.axes() + + plot_data(directory + filename, axs, 'gamma_rho', 'c', '<') + plot_data(directory + filename, axs, 'gamma_u', 'r', 'o') + plot_data(directory + filename, axs, 'gamma_exner', 'b', 's') + plot_data(directory + filename, axs, 'gamma_theta', 'g', '^') + plot_data(directory + filename, axs, 'gamma_mr', 'black', 'x') + + plt.legend(loc='lower right') + plt.xlabel('Gamma') + plt.ylabel('Relative error') + plt.title('Validity of the tangent linear model') + + #plt.show() + plt.savefig(directory + filename + "convergence_plot.png") + + +if __name__ == "__main__": + + DATA_DIRECTORY = os.getcwd()+'/' + DATA_FILENAME = 'outfile' + CONFIG = os.getenv('CONFIG') + print(CONFIG) + + make_plot(DATA_DIRECTORY, DATA_FILENAME) diff --git a/science/linear/plot_convergence/plot_convergence.sh b/science/linear/plot_convergence/plot_convergence.sh new file mode 100755 index 00000000..b0ba236f --- /dev/null +++ b/science/linear/plot_convergence/plot_convergence.sh @@ -0,0 +1,112 @@ +############################################################################## +# (c) Crown copyright 2022 Met Office. All rights reserved. +# The file LICENCE, distributed with this code, contains details of the terms +# under which the code may be used. +############################################################################## + +#---------------------------------------------------------------------- +# Plots the convergence rate of the tangent linear model. +#---------------------------------------------------------------------- + +# INSTRUCTIONS TO RUN +# 1. Specify the config_list +# 2. Run using . plot_convergence.sh, from the plot_convergence directory + +# SCIENCE DETAILS +# The relative linearisation error is +# E = || N(x+ gamma x') - N(x) - L(x) gamma x' || / || L(x) gamma x' || +# where N=nonlinear model, L=linear model, x=linearisation state +# x'=perturbation, gamma=scalar. +# From the Taylor series expansion, E(gamma) = O(gamma) i.e. of the order gamma +# So the relative error should be a linear function of gamma + +# SCRIPT STEPS +# 1. Build the model +# 2. Produce the data: The integration test tl_test_timesteps is extended by +# running over 10 values of gamma, rather than 2 values of gamma. +# 3. Plot the data: The data is plotted for each prognostic variable. + +# EXTENSION +# The plot_configuration.nml can also be extended to other configurations e.g +# * increase the number of timesteps (timesteps_end) +# * increase the number of timesteps between updating the linearisation state +# (update_ls_frequency) + +#-------------------------------------------------------------------------- + +######################## Functions ########################################## +build(){ + echo $CONFIG " Building the executable" + + # Integration tests executable name + exe=$Linear_dir/test/$CONFIG/$CONFIG + + # Build the integration tests, unless that has already been completed + if [ -f $exe ] ; then + echo "Do not need to build the executable as $exe exists" + else + echo "$exe does not exist, so now building the executable" + cd ../../../build + python3 local_build.py -a linear -t integration-tests + + if [$? -ne 0 ]; then + echo "Error building the executable" + return + fi + fi +} + +integration_test(){ + # Setup the configuration - to test with 10 values of gamma + echo $CONFIG " Setting up the configuration" + cd $Linear_dir/integration-test/$CONFIG/resources/ + cp ${CONFIG}_configuration.nml plot_configuration.nml + sed -i 's/number_gamma_values=2/number_gamma_values=10/g' plot_configuration.nml + if [ $? -ne 0 ]; then + echo "Error in creating plot_configuration.nml" + return + else + echo "plot_configuration.nml created in " $PWD + fi + + # Run the tl_test_timesteps integration test + echo $CONFIG " Running the integration test" + cd .. + echo $PWD + $exe resources/plot_configuration.nml test_timesteps > outfile + if [ $? -ne 0 ]; then + echo "Error in creating outfile data" + return + else + echo "Data created successfully" + fi +} + +plot(){ + # Plot the data, together with the expected gradient + echo $CONFIG " Plotting the data" + cd $Linear_dir/integration-test/$CONFIG + python $Linear_dir/plot_convergence/plot_convergence.py +} + +#################### MAIN PROGRAM ################################################# + +# Modules +module purge +module use /home/users/lfricadmin/lmod +module load lfric/vn3.0 +module load scitools + +config_list=(nwp_gal9 semi_implicit runge_kutta) + +# Define directories using the current working directory +export Linear_dir="$(dirname $PWD)" + +for configuration in "${config_list[@]}"; do + echo $configuration + export CONFIG="$configuration" + + build + integration_test + plot +done From 43039e715a3ba8654847a330093400f5f0b0b45d Mon Sep 17 00:00:00 2001 From: Christine Johnson <74597224+cjohnson-pi@users.noreply.github.com> Date: Wed, 21 Jan 2026 11:08:50 +0000 Subject: [PATCH 2/8] integration tests pass for local build --- .../tl_test/tl_test_timesteps_random_alg_mod.x90 | 2 +- science/linear/plot_convergence/plot_convergence.sh | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/science/linear/integration-test/tl_test/tl_test_timesteps_random_alg_mod.x90 b/science/linear/integration-test/tl_test/tl_test_timesteps_random_alg_mod.x90 index 9b8ed013..1d77a9ac 100644 --- a/science/linear/integration-test/tl_test/tl_test_timesteps_random_alg_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_timesteps_random_alg_mod.x90 @@ -224,7 +224,7 @@ module tl_test_timesteps_random_alg_mod gamma_exner = 6.5e-1_r_def if ( use_moisture ) then - gamma_mr = 1.e-2_r_def + gamma_mr = 1.e-3_r_def else gamma_mr = 0.0_r_def end if diff --git a/science/linear/plot_convergence/plot_convergence.sh b/science/linear/plot_convergence/plot_convergence.sh index b0ba236f..1a78a97f 100755 --- a/science/linear/plot_convergence/plot_convergence.sh +++ b/science/linear/plot_convergence/plot_convergence.sh @@ -47,7 +47,7 @@ build(){ else echo "$exe does not exist, so now building the executable" cd ../../../build - python3 local_build.py -a linear -t integration-tests + python3 local_build.py linear -t integration-tests if [$? -ne 0 ]; then echo "Error building the executable" From f6ff3581ab9bee3d84dc82b9b5bc45370b2396ad Mon Sep 17 00:00:00 2001 From: Christine Johnson <74597224+cjohnson-pi@users.noreply.github.com> Date: Wed, 21 Jan 2026 11:19:00 +0000 Subject: [PATCH 3/8] CLA --- CONTRIBUTORS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index d0f7ae14..2b624fda 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -4,3 +4,4 @@ | ----------- | --------- | ----------- | ---- | | james-bruten-mo | James Bruten | Met Office | 2025-12-09 | | jennyhickson | Jenny Hickson | Met Office | 2025-12-10 | +| cjohnson-pi | Christine Johnson | Met Office | 2026-01-21 | From 650bdf059d3fb5baeb773e292b9848fbf8599d46 Mon Sep 17 00:00:00 2001 From: Christine Johnson <74597224+cjohnson-pi@users.noreply.github.com> Date: Wed, 21 Jan 2026 11:53:27 +0000 Subject: [PATCH 4/8] style --- science/gungho/source/driver/gungho_step_mod.x90 | 2 +- .../integration-test/tl_test/tl_test_semi_imp_alg_mod.x90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/science/gungho/source/driver/gungho_step_mod.x90 b/science/gungho/source/driver/gungho_step_mod.x90 index db29183f..8c302dab 100644 --- a/science/gungho/source/driver/gungho_step_mod.x90 +++ b/science/gungho/source/driver/gungho_step_mod.x90 @@ -123,7 +123,7 @@ module gungho_step_mod call log_event( log_scratch_space, LOG_LEVEL_INFO ) temp_corr_io_value => get_io_value( modeldb%values, 'temperature_correction_io_value') - + if ( encorr_usage /= encorr_usage_none ) then call modeldb%values%get_value( 'total_dry_mass', total_dry_mass ) call modeldb%values%get_value( 'total_energy', total_energy ) diff --git a/science/linear/integration-test/tl_test/tl_test_semi_imp_alg_mod.x90 b/science/linear/integration-test/tl_test/tl_test_semi_imp_alg_mod.x90 index 6a58e383..c47a6356 100644 --- a/science/linear/integration-test/tl_test/tl_test_semi_imp_alg_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_semi_imp_alg_mod.x90 @@ -242,9 +242,9 @@ module tl_test_semi_imp_alg_mod do i = 1, nummr call invoke( & setval_X( r_mr(i), theta ), & - inc_a_times_X( 0.000001, r_mr(i) ), & + inc_a_times_X( 0.000001_r_def, r_mr(i) ), & setval_X( ls_mr(i), ls_theta ), & - inc_a_times_X( 0.000001, ls_mr(i) )) + inc_a_times_X( 0.000001_r_def, ls_mr(i) )) end do ! Set the prognostic fields to be the linearisation state From 50bc604f16e35da3427b60029beda8533c669019 Mon Sep 17 00:00:00 2001 From: Christine Johnson <74597224+cjohnson-pi@users.noreply.github.com> Date: Wed, 21 Jan 2026 12:11:06 +0000 Subject: [PATCH 5/8] style --- .../integration-test/tl_test/tl_test_timesteps_alg_mod.x90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/science/linear/integration-test/tl_test/tl_test_timesteps_alg_mod.x90 b/science/linear/integration-test/tl_test/tl_test_timesteps_alg_mod.x90 index 7e5552de..686e0cad 100644 --- a/science/linear/integration-test/tl_test/tl_test_timesteps_alg_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_timesteps_alg_mod.x90 @@ -226,9 +226,9 @@ module tl_test_timesteps_alg_mod do i = 1, nummr call invoke( & setval_X( r_mr(i), theta ), & - inc_a_times_X( 0.000001, r_mr(i) ), & + inc_a_times_X( 0.000001_r_def, r_mr(i) ), & setval_X( ls_mr(i), ls_theta ), & - inc_a_times_X( 0.000001, ls_mr(i) )) + inc_a_times_X( 0.000001_r_def, ls_mr(i) )) end do call moist_dyn_factors_alg( ls_moist_dyn, ls_mr ) From f4bfe9e734f1e9e5782ce670eec0a667e356eca0 Mon Sep 17 00:00:00 2001 From: Christine Johnson <74597224+cjohnson-pi@users.noreply.github.com> Date: Thu, 22 Jan 2026 13:41:14 +0000 Subject: [PATCH 6/8] For gungho_model hot jupiter cases --- science/gungho/source/driver/gungho_step_mod.x90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/science/gungho/source/driver/gungho_step_mod.x90 b/science/gungho/source/driver/gungho_step_mod.x90 index 8c302dab..2354c3a7 100644 --- a/science/gungho/source/driver/gungho_step_mod.x90 +++ b/science/gungho/source/driver/gungho_step_mod.x90 @@ -124,7 +124,7 @@ module gungho_step_mod temp_corr_io_value => get_io_value( modeldb%values, 'temperature_correction_io_value') - if ( encorr_usage /= encorr_usage_none ) then + if ( encorr_usage /= encorr_usage_none .or. write_conservation_diag ) then call modeldb%values%get_value( 'total_dry_mass', total_dry_mass ) call modeldb%values%get_value( 'total_energy', total_energy ) call modeldb%values%get_value( 'total_energy_previous', & From 433b032069bbf5e223579045ba1f6b5c066d2e2b Mon Sep 17 00:00:00 2001 From: Christine Johnson <74597224+cjohnson-pi@users.noreply.github.com> Date: Wed, 28 Jan 2026 14:59:08 +0000 Subject: [PATCH 7/8] science review --- science/gungho/source/driver/gungho_step_mod.x90 | 3 +-- .../tl_test/tl_test_convergence_rate_check.f90 | 14 ++++++++++++-- .../tl_test/tl_test_rhs_alg_mod.x90 | 2 +- .../tl_test/tl_test_semi_imp_alg_mod.x90 | 6 ++++++ .../tl_test/tl_test_timesteps_random_alg_mod.x90 | 8 +++++--- .../linear/plot_convergence/plot_convergence.py | 3 +++ .../linear/plot_convergence/plot_convergence.sh | 8 ++++++-- 7 files changed, 34 insertions(+), 10 deletions(-) diff --git a/science/gungho/source/driver/gungho_step_mod.x90 b/science/gungho/source/driver/gungho_step_mod.x90 index 2354c3a7..82b47579 100644 --- a/science/gungho/source/driver/gungho_step_mod.x90 +++ b/science/gungho/source/driver/gungho_step_mod.x90 @@ -122,9 +122,8 @@ module gungho_step_mod '(A,I0)' ) 'Start of timestep ', model_clock%get_step() call log_event( log_scratch_space, LOG_LEVEL_INFO ) - temp_corr_io_value => get_io_value( modeldb%values, 'temperature_correction_io_value') - if ( encorr_usage /= encorr_usage_none .or. write_conservation_diag ) then + temp_corr_io_value => get_io_value( modeldb%values, 'temperature_correction_io_value') call modeldb%values%get_value( 'total_dry_mass', total_dry_mass ) call modeldb%values%get_value( 'total_energy', total_energy ) call modeldb%values%get_value( 'total_energy_previous', & diff --git a/science/linear/integration-test/tl_test/tl_test_convergence_rate_check.f90 b/science/linear/integration-test/tl_test/tl_test_convergence_rate_check.f90 index 20dee410..b090ba7f 100644 --- a/science/linear/integration-test/tl_test/tl_test_convergence_rate_check.f90 +++ b/science/linear/integration-test/tl_test/tl_test_convergence_rate_check.f90 @@ -20,13 +20,13 @@ module tl_test_convergence_rate_check contains - !> @brief Test the convergence rate. + !> @brief Test the convergence rate, with application of square root. !> @details If the convergence rate is not close to 4, within the !! specified tolerance, set the pass string to FAIL. !> @param[in] name Variable or test name !> @param[in] norm Current norm !> @param[in] norm_prev Previous norm - !> @param[inout] pass_str Pass string (either PASS or FAIL) + !> @param[in,out] pass_str Pass string (either PASS or FAIL) !> @param[in] tol Tolerance subroutine convergence_pass_string( name, norm, norm_prev, pass_str, tol ) implicit none @@ -38,6 +38,16 @@ subroutine convergence_pass_string( name, norm, norm_prev, pass_str, tol ) real(r_def) :: conv_rate + ! Let the error between the nonlinear (N) difference and the linear (L) be + ! E(g) = || N(x + g dx) - N(x) - g Ldx || = O(g^2) + ! where g is a scalar, x and dx are vectors, O is the order + ! and || . || = (x^T x)^1/2 is the L2 norm. + ! + ! Then the ratio + ! E(2g) / E(g) = O(4 g^2) / O(g^2) = 4 + ! i.e. we need to check whether the convergence rate is close to 4. + + ! The norms have not had a square root applied yet. conv_rate = sqrt(norm_prev / norm) if ( abs( conv_rate - 4.0_r_def ) >= tol ) then diff --git a/science/linear/integration-test/tl_test/tl_test_rhs_alg_mod.x90 b/science/linear/integration-test/tl_test/tl_test_rhs_alg_mod.x90 index 794764e8..6207afcd 100644 --- a/science/linear/integration-test/tl_test/tl_test_rhs_alg_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_rhs_alg_mod.x90 @@ -28,7 +28,7 @@ module tl_test_rhs_alg_mod use rhs_alg_mod, only: rhs_alg use tl_rhs_alg_mod, only: tl_rhs_alg use log_mod, only: log_event, LOG_LEVEL_INFO, & - LOG_LEVEL_ERROR + LOG_LEVEL_ERROR use moist_dyn_factors_alg_mod, only: moist_dyn_factors_alg use tl_moist_dyn_factors_alg_mod, only: tl_moist_dyn_factors_alg use tl_test_convergence_rate_check, only: array_convergence_rate_check diff --git a/science/linear/integration-test/tl_test/tl_test_semi_imp_alg_mod.x90 b/science/linear/integration-test/tl_test/tl_test_semi_imp_alg_mod.x90 index c47a6356..30bf004e 100644 --- a/science/linear/integration-test/tl_test/tl_test_semi_imp_alg_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_semi_imp_alg_mod.x90 @@ -239,6 +239,12 @@ module tl_test_semi_imp_alg_mod setval_X(r_rho, rho ) ) ! Overwrite the moisture fields with potential temperature + ! + ! Note: The perturbation fields are read in from file. However the + ! moisture perturbation fields have values close to zero (especially + ! at upper levels) and this gives an inaccurate linearisation test. Using + ! potential temperature is an easy fix to create a realistic 'random' + ! field using real data that has mostly non-zero values. do i = 1, nummr call invoke( & setval_X( r_mr(i), theta ), & diff --git a/science/linear/integration-test/tl_test/tl_test_timesteps_random_alg_mod.x90 b/science/linear/integration-test/tl_test/tl_test_timesteps_random_alg_mod.x90 index 1d77a9ac..6ea6ff29 100644 --- a/science/linear/integration-test/tl_test/tl_test_timesteps_random_alg_mod.x90 +++ b/science/linear/integration-test/tl_test/tl_test_timesteps_random_alg_mod.x90 @@ -292,7 +292,7 @@ module tl_test_timesteps_random_alg_mod ! Clean up solvers and timestep method left over from running ! the model to set up model state and linear state call final_si_operators() - call semi_implicit_solver_alg_final + call semi_implicit_solver_alg_final() call finalise_model(modeldb) call initialise_model( mesh, & @@ -312,7 +312,7 @@ module tl_test_timesteps_random_alg_mod do i = 1, nummr call invoke( setval_X( n1_mr(i), mr(i) ) ) - enddo + end do !------------------------------------------------------------------------ ! Linear model run @@ -323,7 +323,7 @@ module tl_test_timesteps_random_alg_mod setval_X( exner, p_exner) ) do i = 1, nummr call invoke( setval_X( mr(i), p_mr(i) ) ) - enddo + end do call tl_moist_dyn_factors_alg( moist_dyn, mr ) ! Set the linearisation state @@ -525,6 +525,8 @@ module tl_test_timesteps_random_alg_mod !> @brief Test the convergence rate. !> @details If the convergence rate is not close to 4, within the !! specified tolerance, then change the pass string to FAIL. + !! This is a special routine for this particular test, to allow + !! for norms that have already had the square root applied. !> @param[in] norm_diff Current norm !> @param[in] norm_diff_prev Previous norm !> @param[inout] pass_str Pass string (either PASS or FAIL) diff --git a/science/linear/plot_convergence/plot_convergence.py b/science/linear/plot_convergence/plot_convergence.py index e5c7412d..6d7c42ea 100644 --- a/science/linear/plot_convergence/plot_convergence.py +++ b/science/linear/plot_convergence/plot_convergence.py @@ -106,7 +106,10 @@ def make_plot(directory, filename): plt.ylabel('Relative error') plt.title('Validity of the tangent linear model') + # To show the plot to the screen, uncommment plt.show() #plt.show() + + # Save the plot to a file plt.savefig(directory + filename + "convergence_plot.png") diff --git a/science/linear/plot_convergence/plot_convergence.sh b/science/linear/plot_convergence/plot_convergence.sh index 1a78a97f..689a3b9d 100755 --- a/science/linear/plot_convergence/plot_convergence.sh +++ b/science/linear/plot_convergence/plot_convergence.sh @@ -46,7 +46,7 @@ build(){ echo "Do not need to build the executable as $exe exists" else echo "$exe does not exist, so now building the executable" - cd ../../../build + cd $Root_dir/build python3 local_build.py linear -t integration-tests if [$? -ne 0 ]; then @@ -71,7 +71,7 @@ integration_test(){ # Run the tl_test_timesteps integration test echo $CONFIG " Running the integration test" - cd .. + cd $Linear_dir/integration-test/$CONFIG echo $PWD $exe resources/plot_configuration.nml test_timesteps > outfile if [ $? -ne 0 ]; then @@ -101,6 +101,10 @@ config_list=(nwp_gal9 semi_implicit runge_kutta) # Define directories using the current working directory export Linear_dir="$(dirname $PWD)" +export Parent_dir="$(dirname $Linear_dir)" +export Root_dir="$(dirname $Parent_dir)" +echo "Linear_dir" $Linear_dir +echo "Root_dir" $Root_dir for configuration in "${config_list[@]}"; do echo $configuration From 814f038b7cee5bb6325d93076444ce927f671e9a Mon Sep 17 00:00:00 2001 From: Christine Johnson <74597224+cjohnson-pi@users.noreply.github.com> Date: Thu, 29 Jan 2026 11:31:54 +0000 Subject: [PATCH 8/8] science review response --- .../plot_convergence/plot_convergence.py | 81 ------------------ .../plot_convergence/plot_convergence.sh | 83 ------------------- .../plot_convergence/plot_convergence.sh | 6 +- 3 files changed, 4 insertions(+), 166 deletions(-) delete mode 100644 applications/linear_model/plot_convergence/plot_convergence.py delete mode 100755 applications/linear_model/plot_convergence/plot_convergence.sh diff --git a/applications/linear_model/plot_convergence/plot_convergence.py b/applications/linear_model/plot_convergence/plot_convergence.py deleted file mode 100644 index 811421c8..00000000 --- a/applications/linear_model/plot_convergence/plot_convergence.py +++ /dev/null @@ -1,81 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -############################################################################## -# (c) Crown copyright 2021 Met Office. All rights reserved. -# The file LICENCE, distributed with this code, contains details of the terms -# under which the code may be used. -############################################################################## -''' -Using the output derived from the test_timesteps integration test with -10 values of gamma, plot a graph of the relative error (linearisation error) -against the size of the perturbation for different prognostic variables -(gamma). -''' -import os -import pandas as pd -import matplotlib.pyplot as plt - - -def plot_data(filename, axes, variable, color, shape): - ''' - Create a data frame with columns gamma (size of the perturbation) - and norm (the relative error). Plot this data with a log-log axis. - ''' - - datafile = open(filename, 'r') - - norm_line = [] - line = datafile.readline - - while line: - line = datafile.readline() - - if variable in line: - split_line = line.replace('norm', '').replace('\n', '').split('=') - norm_line.append([float(split_line[1]), float(split_line[2])]) - - norm_df = pd.DataFrame(norm_line, columns=['gamma', 'norm']) - - datafile.close() - - norm_df.plot.scatter(x='gamma', y='norm', loglog=True, xlim=(10**0, 10**5), - ylim=(10**-5, 10**0), ax=axes, color=color, - marker=shape) - - -def make_plot(directory, filename): - ''' - Plot the data for the different prognostic variables, together with the - expected gradient ( which is linear ), on the same plot. - ''' - - axs = plt.axes() - - plot_data(directory + filename, axs, 'gamma_rho', 'c', '<') - plot_data(directory + filename, axs, 'gamma_u', 'r', 'o') - plot_data(directory + filename, axs, 'gamma_exner', 'b', 's') - plot_data(directory + filename, axs, 'gamma_theta', 'g', '^') - plot_data(directory + filename, axs, 'gamma_total', 'black', 'x') - plot_data(directory + filename, axs, 'gamma_mr', 'm', '*') - - expected_x = [10**0, 10**5] - expected_y = [10**-5, 10**0] - plt.plot(expected_x, expected_y) - - plt.legend(['Expected gradient (linear)', 'density', 'momentum', - 'exner pressure', 'potential temperature', 'total', - 'moisture mixing ratios'], - loc='lower right') - plt.xlabel('Gamma') - plt.ylabel('Relative error') - plt.title('Validity of the tangent linear model') - - plt.show() - - -if __name__ == "__main__": - - DATA_DIRECTORY = os.getcwd()+'/' - DATA_FILENAME = 'outfile' - - make_plot(DATA_DIRECTORY, DATA_FILENAME) diff --git a/applications/linear_model/plot_convergence/plot_convergence.sh b/applications/linear_model/plot_convergence/plot_convergence.sh deleted file mode 100755 index ab6e366a..00000000 --- a/applications/linear_model/plot_convergence/plot_convergence.sh +++ /dev/null @@ -1,83 +0,0 @@ -############################################################################## -# (c) Crown copyright 2022 Met Office. All rights reserved. -# The file LICENCE, distributed with this code, contains details of the terms -# under which the code may be used. -############################################################################## - -#---------------------------------------------------------------------- -# Plots the convergence rate of the tangent linear model. -#---------------------------------------------------------------------- - -# INSTRUCTIONS TO RUN -# 1. Specify CONFIG -# 2. Run using . plot_convergence.sh, from the plot_convergence directory - -# SCIENCE DETAILS -# The relative linearisation error is -# E = || N(x+ gamma x') - N(x) - L(x) gamma x' || / || L(x) gamma x' || -# where N=nonlinear model, L=linear model, x=linearisation state -# x'=perturbation, gamma=scalar. -# From the Taylor series expansion, E(gamma) = O(gamma) i.e. of the order gamma -# So the relative error should be a linear function of gamma - -# SCRIPT STEPS -# 1. Produce the data: The integration test tl_test_timesteps is extended by -# running over 10 values of gamma, rather than 2 values of gamma. -# 2. Plot the data: The data is plotted for each prognostic variable. - -# EXTENSION -# The plot_configuration.nml can also be extended to other configurations e.g -# * increase the number of timesteps (timesteps_end) -# * increase the number of timesteps between updating the linearisation state -# (update_ls_frequency) - -#-------------------------------------------------------------------------- - -# CONFIG can be specified as either runge_kutta or semi_implicit -CONFIG=semi_implicit - -# Define directories using the current working directory -Working_dir=$PWD -Linear_dir="$(dirname "$PWD")" - -# Integration tests executable name -exe=$Linear_dir/test/$CONFIG - -# Build the integration tests, unless that has already been completed -if [ -f $exe ] ; then - echo "Do not need to build the executable as $exe exists" -else - echo "$exe does not exist, so now building the executable" - cd $Linear_dir - make integration-tests - - if [$? -ne 0 ]; then - echo "Error building the executable" - return - fi -fi - -# Setup the configuration - to test with 10 values of gamma -cd $Linear_dir/test/test_files/$CONFIG -cp ${CONFIG}_configuration.nml plot_configuration.nml -sed -i 's/number_gamma_values=2/number_gamma_values=10/g' plot_configuration.nml -if [ $? -ne 0 ]; then - echo "Error in creating plot_configuration.nml" - return -fi - -# Run the tl_test_timesteps integration test -echo "Running the integration test" -../../$CONFIG plot_configuration.nml test_timesteps > outfile -if [ $? -ne 0 ]; then - echo "Error in creating outfile data" - return -else - echo "Data created successfully" -fi - -# Plot the data, together with the expected gradient -echo "Plotting the data" -python $Working_dir/plot_convergence.py - - diff --git a/science/linear/plot_convergence/plot_convergence.sh b/science/linear/plot_convergence/plot_convergence.sh index 689a3b9d..8b60e533 100755 --- a/science/linear/plot_convergence/plot_convergence.sh +++ b/science/linear/plot_convergence/plot_convergence.sh @@ -99,8 +99,10 @@ module load scitools config_list=(nwp_gal9 semi_implicit runge_kutta) -# Define directories using the current working directory -export Linear_dir="$(dirname $PWD)" +# Directory of this script +SCRIPT_DIR="$(dirname "$(realpath "$BASH_SOURCE")")" +# And parent directories +export Linear_dir="$(dirname $SCRIPT_DIR)" export Parent_dir="$(dirname $Linear_dir)" export Root_dir="$(dirname $Parent_dir)" echo "Linear_dir" $Linear_dir