Fortran examples¶
Note
Compilation script below
Single point¶
This program shows how to compute the temperature and the density for a single point using MCM, DTM and UM. Also, it shows how to get the uncertainties.
Filename: examples/fortran/point.f90
1program point
2
3 use m_um
4 use m_dtm
5 use m_mcm
6
7 implicit none
8
9 ! Path where to find the UM netCDF files in folders 2002, 2004, 2008-2009
10 character(*), parameter :: data_um = "/home/dala/swami/swami/data/um/"
11 ! Path where to find the "DTM_2020_F107_Kp.dat" file
12 character(*), parameter :: data_dtm = "/home/dala/swami/swami/data/"
13
14 real(8) :: altitude = 150.0d0 ! km
15 real(8) :: day_of_year = 53.0d0 ! days
16 real(8) :: local_time = 12.0d0 ! hours
17 real(8) :: latitude = 0d0 ! degrees
18 real(8) :: longitude = 15d0 ! degrees
19 real(8) :: f107 = 140d0 ! F10.7
20 real(8) :: f107m = 139d0 ! F10.7 average
21 real(8) :: kps(2) = [1d0, 1d0] ! Kp
22
23 real(8) :: temp, dens, std_dens, std_temp, unc
24
25 ! Initialise/load the model
26 call init_mcm(data_um, data_dtm)
27
28 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
29 ! MCM
30
31 ! Get the temperature (K) using the MCM model
32 call get_mcm_temp(temp, altitude, latitude, longitude, local_time, day_of_year, &
33 f107, f107m, kps)
34
35 ! Get the density (g/cm3) using the MCM model
36 call get_mcm_dens(dens, altitude, latitude, longitude, local_time, day_of_year, &
37 f107, f107m, kps)
38
39 print *, temp, dens
40
41 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
42 ! UM (up to 152km)
43
44 ! Get the temperature (K) using the UM model
45 call get_um_temp(temp, altitude, latitude, longitude, local_time, day_of_year, &
46 f107, f107m, kps)
47
48 ! Get the density (g/cm3) using the UM model
49 call get_um_dens(dens, altitude, latitude, longitude, local_time, day_of_year, &
50 f107, f107m, kps)
51
52 ! Standard deviation for density (g/cm3) and temperature (K)
53 call get_um_dens_standard_deviation(std_dens, altitude, latitude, longitude, local_time, &
54 day_of_year, f107, f107m, kps)
55 call get_um_temp_standard_deviation(std_temp, altitude, latitude, longitude, local_time, &
56 day_of_year, f107, f107m, kps)
57
58 print *, temp, dens, std_dens, std_temp
59
60 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
61 ! DTM2020 (above 120 km)
62
63 ! Get the temperature (K) density (g/cm3) using the DTM2020 model
64 call get_dtm2020(dens, temp, altitude, latitude, longitude, local_time, day_of_year, &
65 f107, f107m, kps)
66
67 ! Get the uncertainty (%) of the density
68 call get_dtm2020_dens_uncertainty(unc, altitude, latitude, longitude, local_time, &
69 day_of_year, f107, f107m, kps)
70
71 print *, temp, dens, unc
72
73end program point
Altitude profile¶
This program is useful to get the density and temperature profiles at several points.
Filename: examples/fortran/altitude_profile.f90
In examples/fortran it is included a Python script that can plot the generated .dat files.
1program altitude_profile
2
3 use m_um
4 use m_dtm
5 use m_mcm
6
7 implicit none
8
9 character(*), parameter :: data_um = "/home/dala/swami/swami/data/um/"
10 character(*), parameter :: data_dtm = "/home/dala/swami/swami/data/"
11
12 integer :: i
13
14 real(8) :: altitudes(51) = [(dble(i), i=0, 250, 5)]
15 real(8) :: times(2) = [0, 180]
16 real(8) :: local_times(5) = [(dble(i), i=0, 24, 6)]
17 real(8) :: latitudes(7) = [(dble(i), i=-90, 90, 30)]
18 real(8) :: f107(3) = [70d0, 140d0, 170d0]
19 real(8) :: f107m(3) = [70d0, 140d0, 170d0]
20
21 real(8) :: kps(2) = [1d0, 1d0]
22 integer :: ih, it, ilt, ilat, isc
23 character(len=100) :: fname
24 real(8) :: temp, dens
25
26 ! Initialise/load the model (calls lecdtm)
27 call init_mcm(data_um, data_dtm)
28
29 loop_sc: do isc = 1, size(f107) ! solar cycle
30 loop_t: do it = 1, size(times) ! day of the year
31 loop_lt: do ilt = 1, size(local_times) ! local time
32 loop_lat: do ilat = 1, size(latitudes) ! latitude
33 write (fname, '(A,I0.2,A,I0.2,A,I0.2,A,I0.2,A)') "profile_sc", isc, "_lt", ilt, "_lat", ilat, "_t", it, ".dat"
34 open (unit=100, file="altitude_profile/"//"temp_"//fname)
35 open (unit=101, file="altitude_profile/"//"dens_"//fname)
36 write (*, *) "Opening file ", fname
37
38 ! Save metadata in temperature file
39 write (100, *) "# time : ", times(it)
40 write (100, *) "# local_times: ", local_times(ilt)
41 write (100, *) "# latitudes : ", latitudes(ilat)
42 write (100, *) "# f10.7 : ", f107(isc)
43 write (100, *) "# f10.7m : ", f107m(isc)
44 write (100, *) "# kp : ", kps
45
46 ! Save metadata in density file
47 write (101, *) "# time : ", times(it)
48 write (101, *) "# local_times: ", local_times(ilt)
49 write (101, *) "# latitudes : ", latitudes(ilat)
50 write (101, *) "# f10.7 : ", f107(isc)
51 write (101, *) "# f10.7m : ", f107m(isc)
52 write (101, *) "# kp : ", kps
53 loop_h: do ih = 1, size(altitudes)
54
55 ! MCM
56 call get_mcm_temp(temp, altitudes(ih), latitudes(ilat), 0d0, local_times(ilt), times(it), &
57 f107(isc), f107m(isc), kps)
58 call get_mcm_dens(dens, altitudes(ih), latitudes(ilat), 0d0, local_times(ilt), times(it), &
59 f107(isc), f107m(isc), kps)
60
61 ! save results
62 write (100, '(F15.6,F15.6)') altitudes(ih), temp
63 write (101, '(F15.6,E15.6)') altitudes(ih), dens
64 end do loop_h
65 close (100)
66 close (101)
67 end do loop_lat
68 end do loop_lt
69 end do loop_t
70 end do loop_sc
71
72end program altitude_profile
Map at altitude¶
This program is useful to get the density and temperature maps at different altitudes.
Filename: examples/fortran/map_altitude.f90
In examples/fortran it is included a Python script that can plot the generated .dat files.
1program map_altitude
2 use m_mcm, only: get_mcm_dens, get_mcm_temp, init_mcm
3
4 implicit none
5
6 real(8) :: altitude(3) = [90, 140, 300] ! km
7
8 character(*), parameter :: data_um = "/home/dala/swami/swami/data/um/"
9 character(*), parameter :: data_dtm = "/home/dala/swami/swami/data/"
10 character(100) :: fname
11
12 real(8) :: f107 = 100, f107m = 100, kps(2) = [1, 1]
13 real(8) :: doy = 180, aux
14 integer :: i, j, k
15 real(8) :: lat(7) = [(dble(i), i=-90, 90, 30)]
16 real(8) :: lt(5) = [(dble(i), i=0, 24, 6)]
17
18 real(8) :: temp(size(lat), size(lt))
19 real(8) :: dens(size(lat), size(lt))
20
21 ! Load the model
22 call init_mcm(data_um, data_dtm)
23
24 do k = 1, size(altitude)
25 write (fname, '(A,I0.2,A)') "map_", k, ".dat"
26
27 open (unit=100, file="map_altitude/temp_"//trim(fname))
28 open (unit=101, file="map_altitude/dens_"//trim(fname))
29
30 print *, fname
31 temp(:, :) = 0d0
32 dens(:, :) = 0d0
33
34 ! Save metadata in temperature file
35 write (100, *) "# altitude : ", altitude(k)
36 write (100, *) "# latitudes : ", lat
37 write (100, *) "# local_times : ", lt
38 write (100, *) "# f10.7 : ", f107
39 write (100, *) "# f10.7m : ", f107m
40 write (100, *) "# doy : ", doy
41 write (100, *) "# kp : ", kps
42
43 ! Save metadata in density file
44 write (101, *) "# altitude : ", altitude(k)
45 write (101, *) "# latitudes : ", lat
46 write (101, *) "# local_times : ", lt
47 write (101, *) "# f10.7 : ", f107
48 write (101, *) "# f10.7m : ", f107m
49 write (101, *) "# kp : ", kps
50 write (101, *) "# doy : ", doy
51
52 do i = 1, size(lat) ! loop latitude
53 do j = 1, size(lt) ! loop local time
54 call get_mcm_dens(aux, altitude(k), lat(i), 0d0, lt(j), doy, f107, f107m, kps)
55 dens(i, j) = aux
56 call get_mcm_temp(aux, altitude(k), lat(i), 0d0, lt(j), doy, f107, f107m, kps)
57 temp(i, j) = aux
58 end do
59 write (100, *) temp(i, :)
60 write (101, *) dens(i, :)
61 end do
62 close (100)
63 close (101)
64
65 end do
66
67end program map_altitude
UM Winds¶
This example shows how to get the UM winds at a single point.
Filename: examples/fortran/winds.f90
1program winds
2 use m_um, only: get_um_xwind, get_um_ywind, init_um
3
4 implicit none
5
6 ! Path where to find the UM netCDF files in folders 2002, 2004, 2008-2009
7 character(*), parameter :: data_um = "/home/dala/swami/swami/data/um/"
8 ! Path where to find the "DTM_2020_F107_Kp.dat" file
9 character(*), parameter :: data_dtm = "/home/dala/swami/swami/data/"
10
11 real(8) :: altitude = 60.0d0 ! km
12 real(8) :: day_of_year = 53.0d0 ! days
13 real(8) :: local_time = 12.0d0 ! hours
14 real(8) :: latitude = 0d0 ! degrees
15 real(8) :: longitude = 15d0 ! degrees
16 real(8) :: f107 = 140d0 ! F10.7
17 real(8) :: kps(2) = [1d0, 1d0] ! Kp
18 real(8) :: f107m
19
20 real(8) :: x_wind, y_wind
21
22 ! Init model (load path)
23 call init_um(data_um)
24
25
26 ! Solar cycle 1: f107m < 120
27 print*, "Solar cycle 1: f107m < 120"
28 f107m = 70d0
29 call get_um_xwind(x_wind, altitude, latitude, longitude, local_time, day_of_year, f107, f107m, kps)
30 call get_um_ywind(y_wind, altitude, latitude, longitude, local_time, day_of_year, f107, f107m, kps)
31 print*, "f107m = ", f107m, "x_wind = ", x_wind, "y_wind = ", y_wind
32
33 ! Solar cycle 2: 120 < f107m < 160
34 print*, "Solar cycle 2: 120 < f107m < 160"
35 f107m = 140d0
36 call get_um_xwind(x_wind, altitude, latitude, longitude, local_time, day_of_year, f107, f107m, kps)
37 call get_um_ywind(y_wind, altitude, latitude, longitude, local_time, day_of_year, f107, f107m, kps)
38 print*, "f107m = ", f107m, "x_wind = ", x_wind, "y_wind = ", y_wind
39
40 ! Solar cycle 3: f107m > 160
41 print*, "Solar cycle 3: f107m > 160"
42 f107m = 170d0
43 call get_um_xwind(x_wind, altitude, latitude, longitude, local_time, day_of_year, f107, f107m, kps)
44 call get_um_ywind(y_wind, altitude, latitude, longitude, local_time, day_of_year, f107, f107m, kps)
45 print*, "f107m = ", f107m, "x_wind = ", x_wind, "y_wind = ", y_wind
46
47
48end program winds
Compilation script¶
1#!/usr/bin/env bash
2
3SRC='../../src/libswamif'
4FLAGS='-Wall -pedantic -Warray-bounds -fbacktrace'
5# FLAGS=''
6L_NETCDF=$(nf-config --fflags --flibs)
7M_INTERP="$SRC/m_interp.f90"
8M_UM="$M_INTERP $SRC/m_um.f90"
9M_DTM="$SRC/dtm2020_F107_Kp-subr_MCM.f90 $SRC/dtm2020_sigma_function.f90 $SRC/m_dtm.f90"
10M_MCM="$M_UM $M_DTM $SRC/m_mcm.f90"
11
12echo ">>> Delete old tests"
13rm *.o
14rm *.mod
15rm *.x
16rm *.log
17rm *.test
18rm *.dat
19rm *.png
20
21
22echo ">>> SINGLE POINT"
23gfortran $M_MCM point.f90 -o point.x $FLAGS $L_NETCDF && echo "Compilation point succeded!"
24# ./point.x
25
26
27echo ">>> ALTITUDE PROFILE"
28mkdir -p altitude_profile
29gfortran $M_MCM altitude_profile.f90 -o altitude_profile.x $FLAGS $L_NETCDF && echo "Compilation altitude_profiles succeded!"
30# ./altitude_profile.x
31
32echo ">>> map ALTITUDE"
33mkdir -p map_altitude
34gfortran $M_MCM map_altitude.f90 -o map_altitude.x $FLAGS $L_NETCDF && echo "Compilation map_altitude succeded!"
35# ./map_altitude.x
36
37echo ">>> WINDS"
38gfortran $M_MCM winds.f90 -o winds.x $FLAGS $L_NETCDF && echo "Compilation winds succeded!"
39# ./winds.x