fortran协同工作
例子:TELSEM的使用
Fortran的协同工作办法有:
方法一:将多个.f90文件一块编译
gfortran test1.f90 test2.f90 -o exe
方法二:建立静态链接库.lib或动态链接库.dll。需要被链接的库的rwx权限
可参考: 臭石头和雪球的教程:静态库和动态链接库, fortran-hdf库, fortran-netcdf库
方法一.例子.TELSEM的使用
TELSEM由CNRS的F.Aires和C. Priegent 在《A Tool to Estimate Land-Surface Emissivities at Microwave frequencies (TELSEM) for use in numerical weather prediction》 中提出。 TELSEM基于SSM/I地表发射率数据集(C. Prigent, 1997)的19V&H, 22V, 37V&H, 85V&H等通道的Emissivity进行频率插值,得到月平均的,具有一定空间分辨率的(或固定经纬度)的全球陆面微波发射率, 其适用频段为10 – 700 GHz. (cited from胡继恒)。
paper1: Microwave land surface emissivities estimated from SSM/I observations. J. Geophys. Res. Space Phys. 1997
paper2: Land Surface Microwave Emissivities over the Globe for a Decade. Bull. Am. Meteorol. Soc. 2006: 19-85 GHz
paper3: Evaluation of modeled microwave land surface emissivities with satellite-based estimates. J. Geophys. Res. Atmos. 2015: less than 40GHz (JGR-Atmosphere). https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2014JD021817
TELSEM在RTTOV中被用于计算地表发射率,也可以单独使用。TELSEM官网下载包, 解压:
$ cd /home/hjh/rttov13/emis_data
$ tar -xvf telsem2_mw_atlas.tar.bz2
解压后的文件中,mod_mwatlas_m2.F90提供了使用telsem的接口,test_telsem2.F90提供了示例脚本。
胡继恒师兄的目录 /home/hjh/rttov13/emis_data中已解压完毕。我在使用时,只能修改自己的目录文件/data04/1/mlt/run/TELSEM/TELSEM166.f90
不妨直接联合编译和运行:
# 在目录/data04/1/mlt/run/TELSEM下
$h5fc TELSEM166.f90 /home/hjh/rttov13/emis_data/telsem2/mod_mwatlas_m2.F90 -o test.exe
$./test.exe
附: TELSEM166.f90源代码(根据示例test_telsum2.F90修改,模拟了89VH和166VH的地表反射率)
PROGRAM globe_mlse_telsem
! call TELSEM2 atlas and interpolator to generate global 0.25 degree land surface emissivity on monthly basis
! By Qingyang Liu, Aug, 2023, USTC
! adapted from test_telsem.F90
USE mod_mwatlas_m2
use HDF5
IMPLICIT NONE
INTEGER :: m,n
CHARACTER(len=130) fout
CHARACTER(len=2) mm
INTEGER(jpim) :: error_status
LOGICAL(jplm) :: verbose ! For atlas reading subroutine
INTEGER(jpim) :: verb !=1 for TRUE and 0 for FALSE - for emissivity routines
INTEGER(jpim) :: month ! (1->12)
CHARACTER(LEN=256) :: dir ! directory of emis database
REAL(jprb) :: resol ! horizontal resolution for the user
REAL(jprb) :: lat ! (-90->90)
REAL(jprb) :: lon ! (0->360)
REAL(jprb) :: theta ! (0->60�)
! For individual freq interpolations
REAL(jprb) :: ev, eh, stdv, stdh, covvh
REAL(jprb) :: freq ! (19->85GHz)
INTEGER(jpim) :: i
! For multiple freq interpolations
INTEGER(jpim), PARAMETER :: nchan = 2
REAL(jprb) :: ev2(nchan), eh2(nchan), std(2*nchan,2*nchan)
REAL(jprb) :: freq2(nchan)
! Structure containing atlas data
TYPE(telsem2_atlas_data) :: atlas
! !MLSE output
! integer(HID_T) :: fid,dset_id,space,dcpl
! integer :: error
! integer(hsize_t),parameter :: dimo5(3)=(/4,100,50/) ! ouput: Emiss& Tb
! real(jprb) :: MLSEout(4,100,50)!1440,720)
!--- End of header ----------------------------------
verbose = .TRUE. ! Verbose output for reading subroutine
verb = 0 ! No verbose output for emissivity subroutines
!====================================================
! Read the atlas
!====================================================
dir = '/home/hjh/rttov13/emis_data/' ! liting 20240112
month = 6 !! modify or do loop to simulate other months
! WRITE(0,'(a,i3)') 'Reading atlas for month ',month
CALL rttov_readmw_atlas(TRIM(dir), month, atlas, verbose, error_status)
IF (error_status /= 0) THEN
WRITE(0,'(a)') 'Error reading atlas'
STOP 1
ENDIF
!====================================================
! Calculate emissivities
!====================================================
! Multiple frequencies, spatially averaged, return emissivities only
! write(mm,"%0.2i") month !! to the form '06'
! yue = (/"01","02","03","04","05","06","07","08","09","10","11","12"/)
fout ='TELSEM_89GHz_166GHz_Jun.txt'
resol = 0.25
theta = 52.8 !53.1
freq = 10.65
! freq2(:) = (/10.65, 18.7, 23.8, 36.5, 89.0/) ! For multiple frequency experiment
freq2(:) = (/89.0,166.5/) ! For multiple frequency experiment
i = 0 ! Index of freq in freq2 array
open(15,file= trim(adjustl(fout))) !!final product
lat = -90.+0.125
do m = 1,720
lon = 0.125
do n = 1,1440
CALL emis_interp_int_mult(lat, lon, resol, theta, freq2, nchan, atlas, ev2, eh2, verb = verb)
WRITE(15,100) lat,lon,(ev2(i),i=1,nchan),(eh2(i),i=1,nchan)
lon = lon+0.25
end do
lat = lat+0.25
! write(15,100)lat,lon,(ev2(i),i=1,nchan),(eh2(i),i=1,nchan)
print*,m
! MLSEout(:2,n,m)=ev2(:)
! MLSEout(3:,n,m)=eh2(:)
end do
CLOSE(15)
!====================================================
! Close the atlas
!====================================================
CALL rttov_closemw_atlas(atlas)
100 FORMAT(12(1x,f14.6))
! call h5fcreate_f("/data04/1/mlt/run/TELSEM/mlse.hdf",H5F_ACC_TRUNC_F,fid,error)
! call h5screate_simple_f(3, dimo5, space, error)
! PRINT*,'BEGIN'
! call h5pcreate_f(H5P_DATASET_CREATE_F, dcpl, error)
! PRINT*,'BEGIN2'
! call h5dcreate_f(fid, "MLSE", H5T_IEEE_F32LE, space, dset_id, error, dcpl)
! PRINT*,'BEGIN3'
! call h5dwrite_f(dset_id, H5T_IEEE_F32LE, MLSEout, dimo5, error)
! PRINT*,'BEGIN4'
! CALL h5pclose_f(dcpl , error)
! CALL h5dclose_f(dset_id , error)
! CALL h5sclose_f(space, error)
END PROGRAM globe_mlse_telsem
