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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
program point

    use m_um
    use m_dtm
    use m_mcm

    implicit none

    ! Path where to find the UM netCDF files in folders 2002, 2004, 2008-2009
    character(*), parameter :: data_um = "/home/dala/swami/swami/data/um/"
    ! Path where to find the "DTM_2020_F107_Kp.dat" file
    character(*), parameter :: data_dtm = "/home/dala/swami/swami/data/"

    real(8) :: altitude = 150.0d0   ! km
    real(8) :: day_of_year = 53.0d0    ! days
    real(8) :: local_time = 12.0d0    ! hours
    real(8) :: latitude = 0d0       ! degrees
    real(8) :: longitude = 15d0      ! degrees
    real(8) :: f107 = 140d0     ! F10.7
    real(8) :: f107m = 139d0     ! F10.7 average
    real(8) :: kps(2) = [1d0, 1d0]      ! Kp

    real(8) :: temp, dens, std_dens, std_temp, unc

    ! Initialise/load the model
    call init_mcm(data_um, data_dtm)

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! MCM

    ! Get the temperature (K) using the MCM model
    call get_mcm_temp(temp, altitude, latitude, longitude, local_time, day_of_year, &
                      f107, f107m, kps)

    ! Get the density (g/cm3) using the MCM model
    call get_mcm_dens(dens, altitude, latitude, longitude, local_time, day_of_year, &
                      f107, f107m, kps)

    print *, temp, dens

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! UM  (up to 152km)

    ! Get the temperature (K) using the UM model
    call get_um_temp(temp, altitude, latitude, longitude, local_time, day_of_year, &
                     f107, f107m, kps)

    ! Get the density (g/cm3) using the UM model
    call get_um_dens(dens, altitude, latitude, longitude, local_time, day_of_year, &
                     f107, f107m, kps)

    ! Standard deviation for density (g/cm3) and temperature (K)
    call get_um_dens_standard_deviation(std_dens, altitude, latitude, longitude, local_time, &
                                        day_of_year, f107, f107m, kps)
    call get_um_temp_standard_deviation(std_temp, altitude, latitude, longitude, local_time, &
                                        day_of_year, f107, f107m, kps)

    print *, temp, dens, std_dens, std_temp

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! DTM2020 (above 120 km)

    ! Get the temperature (K) density (g/cm3) using the DTM2020 model
    call get_dtm2020(dens, temp, altitude, latitude, longitude, local_time, day_of_year, &
                     f107, f107m, kps)

    ! Get the uncertainty (%) of the density
    call get_dtm2020_dens_uncertainty(unc, altitude, latitude, longitude, local_time, &
                                      day_of_year, f107, f107m, kps)

    print *, temp, dens, unc

end 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.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
program altitude_profile

    use m_um
    use m_dtm
    use m_mcm

    implicit none

    character(*), parameter :: data_um = "/home/dala/swami/swami/data/um/"
    character(*), parameter :: data_dtm = "/home/dala/swami/swami/data/"

    integer :: i

    real(8) :: altitudes(51) = [(dble(i), i=0, 250, 5)]
    real(8) :: times(2) = [0, 180]
    real(8) :: local_times(5) = [(dble(i), i=0, 24, 6)]
    real(8) :: latitudes(7) = [(dble(i), i=-90, 90, 30)]
    real(8) :: f107(3) = [70d0, 140d0, 170d0]
    real(8) :: f107m(3) = [70d0, 140d0, 170d0]

    real(8) :: kps(2) = [1d0, 1d0]
    integer :: ih, it, ilt, ilat, isc
    character(len=100) :: fname
    real(8) :: temp, dens

    ! Initialise/load the model (calls lecdtm)
    call init_mcm(data_um, data_dtm)

    loop_sc: do isc = 1, size(f107)         ! solar cycle
    loop_t: do it = 1, size(times)           ! day of the year
    loop_lt: do ilt = 1, size(local_times)   ! local time
    loop_lat: do ilat = 1, size(latitudes)   ! latitude
        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"
        open (unit=100, file="altitude_profile/"//"temp_"//fname)
        open (unit=101, file="altitude_profile/"//"dens_"//fname)
        write (*, *) "Opening file ", fname

        ! Save metadata in temperature file
        write (100, *) "# time       : ", times(it)
        write (100, *) "# local_times: ", local_times(ilt)
        write (100, *) "# latitudes  : ", latitudes(ilat)
        write (100, *) "# f10.7      : ", f107(isc)
        write (100, *) "# f10.7m     : ", f107m(isc)
        write (100, *) "# kp         : ", kps

        ! Save metadata in density file
        write (101, *) "# time       : ", times(it)
        write (101, *) "# local_times: ", local_times(ilt)
        write (101, *) "# latitudes  : ", latitudes(ilat)
        write (101, *) "# f10.7      : ", f107(isc)
        write (101, *) "# f10.7m     : ", f107m(isc)
        write (101, *) "# kp         : ", kps
        loop_h: do ih = 1, size(altitudes)

            ! MCM
            call get_mcm_temp(temp, altitudes(ih), latitudes(ilat), 0d0, local_times(ilt), times(it), &
                              f107(isc), f107m(isc), kps)
            call get_mcm_dens(dens, altitudes(ih), latitudes(ilat), 0d0, local_times(ilt), times(it), &
                              f107(isc), f107m(isc), kps)

            ! save results
            write (100, '(F15.6,F15.6)') altitudes(ih), temp
            write (101, '(F15.6,E15.6)') altitudes(ih), dens
        end do loop_h
        close (100)
        close (101)
    end do loop_lat
    end do loop_lt
    end do loop_t
    end do loop_sc

end 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.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
program map_altitude
    use m_mcm, only: get_mcm_dens, get_mcm_temp, init_mcm

    implicit none

    real(8) :: altitude(3) = [90, 140, 300]  ! km

    character(*), parameter :: data_um = "/home/dala/swami/swami/data/um/"
    character(*), parameter :: data_dtm = "/home/dala/swami/swami/data/"
    character(100) :: fname

    real(8) :: f107 = 100, f107m = 100, kps(2) = [1, 1]
    real(8) :: doy = 180, aux
    integer :: i, j, k
    real(8) :: lat(7) = [(dble(i), i=-90, 90, 30)]
    real(8) :: lt(5) = [(dble(i), i=0, 24, 6)]

    real(8) :: temp(size(lat), size(lt))
    real(8) :: dens(size(lat), size(lt))

    ! Load the model
    call init_mcm(data_um, data_dtm)

    do k = 1, size(altitude)
        write (fname, '(A,I0.2,A)') "map_", k, ".dat"

        open (unit=100, file="map_altitude/temp_"//trim(fname))
        open (unit=101, file="map_altitude/dens_"//trim(fname))

        print *, fname
        temp(:, :) = 0d0
        dens(:, :) = 0d0

        ! Save metadata in temperature file
        write (100, *) "# altitude : ", altitude(k)
        write (100, *) "# latitudes : ", lat
        write (100, *) "# local_times : ", lt
        write (100, *) "# f10.7 : ", f107
        write (100, *) "# f10.7m : ", f107m
        write (100, *) "# doy : ", doy
        write (100, *) "# kp : ", kps

        ! Save metadata in density file
        write (101, *) "# altitude : ", altitude(k)
        write (101, *) "# latitudes : ", lat
        write (101, *) "# local_times : ", lt
        write (101, *) "# f10.7 : ", f107
        write (101, *) "# f10.7m : ", f107m
        write (101, *) "# kp : ", kps
        write (101, *) "# doy : ", doy

        do i = 1, size(lat)     ! loop latitude
            do j = 1, size(lt)  ! loop local time
                call get_mcm_dens(aux, altitude(k), lat(i), 0d0, lt(j), doy, f107, f107m, kps)
                dens(i, j) = aux
                call get_mcm_temp(aux, altitude(k), lat(i), 0d0, lt(j), doy, f107, f107m, kps)
                temp(i, j) = aux
            end do
            write (100, *) temp(i, :)
            write (101, *) dens(i, :)
        end do
        close (100)
        close (101)

    end do

end program map_altitude

UM Winds

This example shows how to get the UM winds at a single point.

Filename: examples/fortran/winds.f90

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
program winds
    use m_um, only: get_um_xwind, get_um_ywind, init_um

    implicit none

    ! Path where to find the UM netCDF files in folders 2002, 2004, 2008-2009
    character(*), parameter :: data_um = "/home/dala/swami/swami/data/um/"
    ! Path where to find the "DTM_2020_F107_Kp.dat" file
    character(*), parameter :: data_dtm = "/home/dala/swami/swami/data/"

    real(8) :: altitude = 60.0d0   ! km
    real(8) :: day_of_year = 53.0d0    ! days
    real(8) :: local_time = 12.0d0    ! hours
    real(8) :: latitude = 0d0       ! degrees
    real(8) :: longitude = 15d0      ! degrees
    real(8) :: f107 = 140d0     ! F10.7
    real(8) :: kps(2) = [1d0, 1d0]      ! Kp
    real(8) :: f107m

    real(8) :: x_wind, y_wind

    ! Init model (load path)
    call init_um(data_um)


    ! Solar cycle 1: f107m < 120
    print*, "Solar cycle 1: f107m < 120"
    f107m = 70d0
    call get_um_xwind(x_wind, altitude, latitude, longitude, local_time, day_of_year, f107, f107m, kps)
    call get_um_ywind(y_wind, altitude, latitude, longitude, local_time, day_of_year, f107, f107m, kps)
    print*, "f107m = ", f107m, "x_wind = ", x_wind, "y_wind = ", y_wind
    
    ! Solar cycle 2: 120 < f107m < 160
    print*, "Solar cycle 2: 120 < f107m < 160"
    f107m = 140d0
    call get_um_xwind(x_wind, altitude, latitude, longitude, local_time, day_of_year, f107, f107m, kps)
    call get_um_ywind(y_wind, altitude, latitude, longitude, local_time, day_of_year, f107, f107m, kps)
    print*, "f107m = ", f107m, "x_wind = ", x_wind, "y_wind = ", y_wind
    
    ! Solar cycle 3: f107m > 160
    print*, "Solar cycle 3: f107m > 160"
    f107m = 170d0
    call get_um_xwind(x_wind, altitude, latitude, longitude, local_time, day_of_year, f107, f107m, kps)
    call get_um_ywind(y_wind, altitude, latitude, longitude, local_time, day_of_year, f107, f107m, kps)
    print*, "f107m = ", f107m, "x_wind = ", x_wind, "y_wind = ", y_wind


end program winds

Compilation script

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#!/usr/bin/env bash

SRC='../../src/libswamif'
FLAGS='-Wall -pedantic -Warray-bounds -fbacktrace'
# FLAGS=''
L_NETCDF=$(nf-config --fflags --flibs)
M_INTERP="$SRC/m_interp.f90"
M_UM="$M_INTERP $SRC/m_um.f90"
M_DTM="$SRC/dtm2020_F107_Kp-subr_MCM.f90 $SRC/dtm2020_sigma_function.f90 $SRC/m_dtm.f90"
M_MCM="$M_UM $M_DTM $SRC/m_mcm.f90"

echo ">>> Delete old tests"
rm *.o
rm *.mod
rm *.x
rm *.log
rm *.test
rm *.dat
rm *.png


echo ">>> SINGLE POINT"
gfortran $M_MCM point.f90 -o point.x $FLAGS $L_NETCDF && echo "Compilation point succeded!"
# ./point.x


echo ">>> ALTITUDE PROFILE"
mkdir -p altitude_profile
gfortran $M_MCM altitude_profile.f90 -o altitude_profile.x $FLAGS $L_NETCDF && echo "Compilation altitude_profiles succeded!"
# ./altitude_profile.x

echo ">>> map ALTITUDE"
mkdir -p map_altitude
gfortran $M_MCM map_altitude.f90 -o map_altitude.x $FLAGS $L_NETCDF && echo "Compilation map_altitude succeded!"
# ./map_altitude.x

echo ">>> WINDS"
gfortran $M_MCM winds.f90 -o winds.x $FLAGS $L_NETCDF && echo "Compilation winds succeded!"
# ./winds.x