[1179] | 1 | ! $Id: readaerosol.F90 2346 2015-08-21 15:13:46Z musat $ |
---|
[524] | 2 | ! |
---|
[1179] | 3 | MODULE readaerosol_mod |
---|
| 4 | |
---|
| 5 | REAL, SAVE :: not_valid=-333. |
---|
| 6 | |
---|
| 7 | CONTAINS |
---|
| 8 | |
---|
[1492] | 9 | SUBROUTINE readaerosol(name_aero, type, filename, iyr_in, klev_src, pt_ap, pt_b, pt_out, psurf, load) |
---|
[1179] | 10 | |
---|
| 11 | !**************************************************************************************** |
---|
| 12 | ! This routine will read the aersosol from file. |
---|
[524] | 13 | ! |
---|
[1179] | 14 | ! Read a year data with get_aero_fromfile depending on aer_type : |
---|
| 15 | ! - actuel : read year 1980 |
---|
| 16 | ! - preind : read natural data |
---|
| 17 | ! - scenario : read one or two years and do eventually linare time interpolation |
---|
[1150] | 18 | ! |
---|
[1179] | 19 | ! Return pointer, pt_out, to the year read or result from interpolation |
---|
| 20 | !**************************************************************************************** |
---|
| 21 | USE dimphy |
---|
[2311] | 22 | USE print_control_mod, ONLY: lunout |
---|
[524] | 23 | |
---|
[1179] | 24 | IMPLICIT NONE |
---|
[766] | 25 | |
---|
[1179] | 26 | ! Input arguments |
---|
| 27 | CHARACTER(len=7), INTENT(IN) :: name_aero |
---|
[1492] | 28 | CHARACTER(len=*), INTENT(IN) :: type ! actuel, annuel, scenario or preind |
---|
| 29 | CHARACTER(len=8), INTENT(IN) :: filename |
---|
[1179] | 30 | INTEGER, INTENT(IN) :: iyr_in |
---|
[1150] | 31 | |
---|
[1179] | 32 | ! Output |
---|
| 33 | INTEGER, INTENT(OUT) :: klev_src |
---|
| 34 | REAL, POINTER, DIMENSION(:) :: pt_ap ! Pointer for describing the vertical levels |
---|
| 35 | REAL, POINTER, DIMENSION(:) :: pt_b ! Pointer for describing the vertical levels |
---|
[1270] | 36 | REAL, POINTER, DIMENSION(:,:,:) :: pt_out ! The massvar distributions, DIMENSION(klon, klev_src, 12) |
---|
| 37 | REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf ! Surface pression for 12 months |
---|
| 38 | REAL, DIMENSION(klon,12), INTENT(OUT) :: load ! Aerosol mass load in each column for 12 months |
---|
[1150] | 39 | |
---|
[1179] | 40 | ! Local variables |
---|
[1270] | 41 | CHARACTER(len=4) :: cyear |
---|
| 42 | REAL, POINTER, DIMENSION(:,:,:) :: pt_2 |
---|
| 43 | REAL, DIMENSION(klon,12) :: psurf2, load2 |
---|
| 44 | REAL :: p0 ! Reference pressure |
---|
| 45 | INTEGER :: iyr1, iyr2, klev_src2 |
---|
| 46 | INTEGER :: it, k, i |
---|
| 47 | LOGICAL, PARAMETER :: lonlyone=.FALSE. |
---|
[1150] | 48 | |
---|
[1179] | 49 | !**************************************************************************************** |
---|
| 50 | ! Read data depending on aer_type |
---|
| 51 | ! |
---|
| 52 | !**************************************************************************************** |
---|
[524] | 53 | |
---|
[1179] | 54 | IF (type == 'actuel') THEN |
---|
| 55 | ! Read and return data for year 1980 |
---|
| 56 | !**************************************************************************************** |
---|
| 57 | cyear='1980' |
---|
[1270] | 58 | ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month |
---|
| 59 | ! pt_out has dimensions (klon, klev_src, 12) |
---|
[1492] | 60 | CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load) |
---|
[1179] | 61 | |
---|
[640] | 62 | |
---|
[1179] | 63 | ELSE IF (type == 'preind') THEN |
---|
| 64 | ! Read and return data from file with suffix .nat |
---|
| 65 | !**************************************************************************************** |
---|
| 66 | cyear='.nat' |
---|
[1270] | 67 | ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month |
---|
| 68 | ! pt_out has dimensions (klon, klev_src, 12) |
---|
[1492] | 69 | CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load) |
---|
[1179] | 70 | |
---|
[1321] | 71 | ELSE IF (type == 'annuel') THEN |
---|
| 72 | ! Read and return data from scenario annual files |
---|
| 73 | !**************************************************************************************** |
---|
| 74 | WRITE(cyear,'(I4)') iyr_in |
---|
| 75 | WRITE(lunout,*) 'get_aero 3 iyr_in=', iyr_in,' ',cyear |
---|
| 76 | ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month |
---|
| 77 | ! pt_out has dimensions (klon, klev_src, 12) |
---|
[1492] | 78 | CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load) |
---|
[1321] | 79 | |
---|
[1179] | 80 | ELSE IF (type == 'scenario') THEN |
---|
| 81 | ! Read data depending on actual year and interpolate if necessary |
---|
| 82 | !**************************************************************************************** |
---|
| 83 | IF (iyr_in .LT. 1850) THEN |
---|
| 84 | cyear='.nat' |
---|
| 85 | WRITE(lunout,*) 'get_aero 1 iyr_in=', iyr_in,' ',cyear |
---|
[1270] | 86 | ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month |
---|
| 87 | ! pt_out has dimensions (klon, klev_src, 12) |
---|
[1492] | 88 | CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load) |
---|
[1179] | 89 | |
---|
| 90 | ELSE IF (iyr_in .GE. 2100) THEN |
---|
| 91 | cyear='2100' |
---|
| 92 | WRITE(lunout,*) 'get_aero 2 iyr_in=', iyr_in,' ',cyear |
---|
[1270] | 93 | ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month |
---|
| 94 | ! pt_out has dimensions (klon, klev_src, 12) |
---|
[1492] | 95 | CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load) |
---|
[1179] | 96 | |
---|
| 97 | ELSE |
---|
| 98 | ! Read data from 2 decades and interpolate to actual year |
---|
| 99 | ! a) from actual 10-yr-period |
---|
| 100 | IF (iyr_in.LT.1900) THEN |
---|
| 101 | iyr1 = 1850 |
---|
| 102 | iyr2 = 1900 |
---|
| 103 | ELSE IF (iyr_in.GE.1900.AND.iyr_in.LT.1920) THEN |
---|
| 104 | iyr1 = 1900 |
---|
| 105 | iyr2 = 1920 |
---|
| 106 | ELSE |
---|
| 107 | iyr1 = INT(iyr_in/10)*10 |
---|
| 108 | iyr2 = INT(1+iyr_in/10)*10 |
---|
| 109 | ENDIF |
---|
| 110 | |
---|
| 111 | WRITE(cyear,'(I4)') iyr1 |
---|
| 112 | WRITE(lunout,*) 'get_aero 3 iyr_in=', iyr_in,' ',cyear |
---|
[1270] | 113 | ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month |
---|
| 114 | ! pt_out has dimensions (klon, klev_src, 12) |
---|
[1492] | 115 | CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load) |
---|
[1179] | 116 | |
---|
| 117 | ! If to read two decades: |
---|
| 118 | IF (.NOT.lonlyone) THEN |
---|
| 119 | |
---|
| 120 | ! b) from the next following one |
---|
| 121 | WRITE(cyear,'(I4)') iyr2 |
---|
| 122 | WRITE(lunout,*) 'get_aero 4 iyr_in=', iyr_in,' ',cyear |
---|
| 123 | |
---|
| 124 | NULLIFY(pt_2) |
---|
[1270] | 125 | ! get_aero_fromfile returns pt_2 allocated and initialized with data for 12 month |
---|
| 126 | ! pt_2 has dimensions (klon, klev_src, 12) |
---|
[1492] | 127 | CALL get_aero_fromfile(name_aero, cyear, filename, klev_src2, pt_ap, pt_b, p0, pt_2, psurf2, load2) |
---|
[1179] | 128 | ! Test for same number of vertical levels |
---|
| 129 | IF (klev_src /= klev_src2) THEN |
---|
| 130 | WRITE(lunout,*) 'Two aerosols files with different number of vertical levels is not allowded' |
---|
[2311] | 131 | CALL abort_physic('readaersosol','Error in number of vertical levels',1) |
---|
[1179] | 132 | END IF |
---|
| 133 | |
---|
| 134 | ! Linare interpolate to the actual year: |
---|
[1270] | 135 | DO it=1,12 |
---|
[1179] | 136 | DO k=1,klev_src |
---|
| 137 | DO i = 1, klon |
---|
| 138 | pt_out(i,k,it) = & |
---|
[1403] | 139 | pt_out(i,k,it) - REAL(iyr_in-iyr1)/REAL(iyr2-iyr1) * & |
---|
[1179] | 140 | (pt_out(i,k,it) - pt_2(i,k,it)) |
---|
| 141 | END DO |
---|
| 142 | END DO |
---|
[1143] | 143 | |
---|
[1179] | 144 | DO i = 1, klon |
---|
[1270] | 145 | psurf(i,it) = & |
---|
[1403] | 146 | psurf(i,it) - REAL(iyr_in-iyr1)/REAL(iyr2-iyr1) * & |
---|
[1270] | 147 | (psurf(i,it) - psurf2(i,it)) |
---|
[640] | 148 | |
---|
[1270] | 149 | load(i,it) = & |
---|
[1403] | 150 | load(i,it) - REAL(iyr_in-iyr1)/REAL(iyr2-iyr1) * & |
---|
[1270] | 151 | (load(i,it) - load2(i,it)) |
---|
[1179] | 152 | END DO |
---|
| 153 | END DO |
---|
[524] | 154 | |
---|
[1179] | 155 | ! Deallocate pt_2 no more needed |
---|
| 156 | DEALLOCATE(pt_2) |
---|
| 157 | |
---|
| 158 | END IF ! lonlyone |
---|
| 159 | END IF ! iyr_in .LT. 1850 |
---|
[524] | 160 | |
---|
[1179] | 161 | ELSE |
---|
[1492] | 162 | WRITE(lunout,*)'This option is not implemented : aer_type = ', type,' name_aero=',name_aero |
---|
[2311] | 163 | CALL abort_physic('readaerosol','Error : aer_type parameter not accepted',1) |
---|
[1179] | 164 | END IF ! type |
---|
[640] | 165 | |
---|
[1151] | 166 | |
---|
[1179] | 167 | END SUBROUTINE readaerosol |
---|
[1143] | 168 | |
---|
| 169 | |
---|
[1492] | 170 | SUBROUTINE get_aero_fromfile(varname, cyr, filename, klev_src, pt_ap, pt_b, p0, pt_year, psurf_out, load_out) |
---|
[1179] | 171 | !**************************************************************************************** |
---|
[1270] | 172 | ! Read 12 month aerosol from file and distribute to local process on physical grid. |
---|
[1179] | 173 | ! Vertical levels, klev_src, may differ from model levels if new file format. |
---|
| 174 | ! |
---|
| 175 | ! For mpi_root and master thread : |
---|
| 176 | ! 1) Open file |
---|
| 177 | ! 2) Find vertical dimension klev_src |
---|
| 178 | ! 3) Read field month by month |
---|
| 179 | ! 4) Close file |
---|
[2346] | 180 | ! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo) |
---|
[1179] | 181 | ! - Also the levels and the latitudes have to be inversed |
---|
| 182 | ! |
---|
| 183 | ! For all processes and threads : |
---|
| 184 | ! 6) Scatter global field(klon_glo) to local process domain(klon) |
---|
| 185 | ! 7) Test for negative values |
---|
| 186 | !**************************************************************************************** |
---|
[1150] | 187 | |
---|
[1179] | 188 | USE netcdf |
---|
| 189 | USE dimphy |
---|
[2346] | 190 | USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, & |
---|
| 191 | grid2Dto1D_glo |
---|
[1179] | 192 | USE mod_phys_lmdz_para |
---|
[1183] | 193 | USE iophy, ONLY : io_lon, io_lat |
---|
[2311] | 194 | USE print_control_mod, ONLY: lunout |
---|
[1179] | 195 | |
---|
| 196 | IMPLICIT NONE |
---|
[524] | 197 | |
---|
[1179] | 198 | ! Input argumets |
---|
| 199 | CHARACTER(len=7), INTENT(IN) :: varname |
---|
| 200 | CHARACTER(len=4), INTENT(IN) :: cyr |
---|
[1492] | 201 | CHARACTER(len=8), INTENT(IN) :: filename |
---|
[1150] | 202 | |
---|
[1179] | 203 | ! Output arguments |
---|
| 204 | INTEGER, INTENT(OUT) :: klev_src ! Number of vertical levels in file |
---|
| 205 | REAL, POINTER, DIMENSION(:) :: pt_ap ! Pointer for describing the vertical levels |
---|
| 206 | REAL, POINTER, DIMENSION(:) :: pt_b ! Pointer for describing the vertical levels |
---|
| 207 | REAL :: p0 ! Reference pressure value |
---|
[1270] | 208 | REAL, POINTER, DIMENSION(:,:,:) :: pt_year ! Pointer-variabale from file, 12 month, grid : klon,klev_src |
---|
| 209 | REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf_out ! Surface pression for 12 months |
---|
| 210 | REAL, DIMENSION(klon,12), INTENT(OUT) :: load_out ! Aerosol mass load in each column |
---|
[1319] | 211 | INTEGER :: nbr_tsteps ! number of month in file read |
---|
[1150] | 212 | |
---|
[1179] | 213 | ! Local variables |
---|
| 214 | CHARACTER(len=30) :: fname |
---|
| 215 | CHARACTER(len=30) :: cvar |
---|
| 216 | INTEGER :: ncid, dimid, varid |
---|
| 217 | INTEGER :: imth, i, j, k, ierr |
---|
| 218 | REAL :: npole, spole |
---|
| 219 | REAL, ALLOCATABLE, DIMENSION(:,:,:) :: varmth |
---|
[1270] | 220 | REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: varyear ! Global variable read from file, 12 month |
---|
| 221 | REAL, ALLOCATABLE, DIMENSION(:,:,:) :: varyear_glo1D !(klon_glo, klev_src, 12) |
---|
[1179] | 222 | REAL, ALLOCATABLE, DIMENSION(:) :: varktmp |
---|
[524] | 223 | |
---|
[2346] | 224 | REAL, DIMENSION(nbp_lon,nbp_lat,12) :: psurf_glo2D ! Surface pression for 12 months on dynamics global grid |
---|
[1270] | 225 | REAL, DIMENSION(klon_glo,12) :: psurf_glo1D ! -"- on physical global grid |
---|
[2346] | 226 | REAL, DIMENSION(nbp_lon,nbp_lat,12) :: load_glo2D ! Load for 12 months on dynamics global grid |
---|
[1270] | 227 | REAL, DIMENSION(klon_glo,12) :: load_glo1D ! -"- on physical global grid |
---|
[2346] | 228 | REAL, DIMENSION(nbp_lon,nbp_lat) :: vartmp |
---|
| 229 | REAL, DIMENSION(nbp_lon) :: lon_src ! longitudes in file |
---|
| 230 | REAL, DIMENSION(nbp_lat) :: lat_src, lat_src_inv ! latitudes in file |
---|
[1183] | 231 | LOGICAL :: new_file ! true if new file format detected |
---|
| 232 | LOGICAL :: invert_lat ! true if the field has to be inverted for latitudes |
---|
[766] | 233 | |
---|
[524] | 234 | |
---|
[1179] | 235 | ! Deallocate pointers |
---|
| 236 | IF (ASSOCIATED(pt_ap)) DEALLOCATE(pt_ap) |
---|
| 237 | IF (ASSOCIATED(pt_b)) DEALLOCATE(pt_b) |
---|
[766] | 238 | |
---|
[1249] | 239 | IF (is_mpi_root .AND. is_omp_root) THEN |
---|
[524] | 240 | |
---|
[1179] | 241 | ! 1) Open file |
---|
| 242 | !**************************************************************************************** |
---|
[1492] | 243 | ! Add suffix to filename |
---|
| 244 | fname = trim(filename)//cyr//'.nc' |
---|
[1150] | 245 | |
---|
[1492] | 246 | WRITE(lunout,*) 'reading variable ',TRIM(varname),' in file ', TRIM(fname) |
---|
[1609] | 247 | CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid), "pb open "//trim(varname) ) |
---|
[1150] | 248 | |
---|
[1183] | 249 | ! Test for equal longitudes and latitudes in file and model |
---|
| 250 | !**************************************************************************************** |
---|
| 251 | ! Read and test longitudes |
---|
[1609] | 252 | CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" ) |
---|
| 253 | CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" ) |
---|
[1183] | 254 | |
---|
[1202] | 255 | IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN |
---|
[1183] | 256 | WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname) |
---|
| 257 | WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src |
---|
| 258 | WRITE(lunout,*) 'longitudes in model :', io_lon |
---|
| 259 | |
---|
[2311] | 260 | CALL abort_physic('get_aero_fromfile', 'longitudes are not the same in file and model',1) |
---|
[1183] | 261 | END IF |
---|
| 262 | |
---|
| 263 | ! Read and test latitudes |
---|
[1609] | 264 | CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" ) |
---|
| 265 | CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" ) |
---|
[1183] | 266 | |
---|
| 267 | ! Invert source latitudes |
---|
[2346] | 268 | DO j = 1, nbp_lat |
---|
| 269 | lat_src_inv(j) = lat_src(nbp_lat +1 -j) |
---|
[1183] | 270 | END DO |
---|
| 271 | |
---|
[1202] | 272 | IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN |
---|
[1183] | 273 | ! Latitudes are the same |
---|
| 274 | invert_lat=.FALSE. |
---|
[1202] | 275 | ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN |
---|
[1183] | 276 | ! Inverted source latitudes correspond to model latitudes |
---|
| 277 | WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname) |
---|
| 278 | invert_lat=.TRUE. |
---|
| 279 | ELSE |
---|
| 280 | WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname) |
---|
| 281 | WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src |
---|
| 282 | WRITE(lunout,*) 'latitudes in model :', io_lat |
---|
[2311] | 283 | CALL abort_physic('get_aero_fromfile', 'latitudes do not correspond between file and model',1) |
---|
[1183] | 284 | END IF |
---|
| 285 | |
---|
[1319] | 286 | |
---|
[1179] | 287 | ! 2) Check if old or new file is avalabale. |
---|
| 288 | ! New type of file should contain the dimension 'lev' |
---|
| 289 | ! Old type of file should contain the dimension 'PRESNIVS' |
---|
| 290 | !**************************************************************************************** |
---|
| 291 | ierr = nf90_inq_dimid(ncid, 'lev', dimid) |
---|
| 292 | IF (ierr /= NF90_NOERR) THEN |
---|
| 293 | ! Coordinate axe lev not found. Check for presnivs. |
---|
| 294 | ierr = nf90_inq_dimid(ncid, 'PRESNIVS', dimid) |
---|
| 295 | IF (ierr /= NF90_NOERR) THEN |
---|
| 296 | ! Dimension PRESNIVS not found either |
---|
[2311] | 297 | CALL abort_physic('get_aero_fromfile', 'dimension lev or presnivs not in file',1) |
---|
[1179] | 298 | ELSE |
---|
| 299 | ! Old file found |
---|
| 300 | new_file=.FALSE. |
---|
| 301 | WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will not be done' |
---|
| 302 | END IF |
---|
| 303 | ELSE |
---|
| 304 | ! New file found |
---|
| 305 | new_file=.TRUE. |
---|
| 306 | WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will be done' |
---|
| 307 | END IF |
---|
| 308 | |
---|
| 309 | ! 2) Find vertical dimension klev_src |
---|
| 310 | !**************************************************************************************** |
---|
[1609] | 311 | CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src),"pb inq dim for PRESNIVS or lev" ) |
---|
[1179] | 312 | |
---|
| 313 | ! Allocate variables depending on the number of vertical levels |
---|
[2346] | 314 | ALLOCATE(varmth(nbp_lon,nbp_lat, klev_src), varyear(nbp_lon,nbp_lat, klev_src, 12), stat=ierr) |
---|
[2311] | 315 | IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 1',1) |
---|
[1150] | 316 | |
---|
[1179] | 317 | ALLOCATE(pt_ap(klev_src), pt_b(klev_src), varktmp(klev_src), stat=ierr) |
---|
[2311] | 318 | IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 2',1) |
---|
[1150] | 319 | |
---|
[1179] | 320 | ! 3) Read all variables from file |
---|
| 321 | ! There is 2 options for the file structure : |
---|
| 322 | ! new_file=TRUE : read varyear, ps, pt_ap and pt_b |
---|
| 323 | ! new_file=FALSE : read varyear month by month |
---|
| 324 | !**************************************************************************************** |
---|
[1150] | 325 | |
---|
[1179] | 326 | IF (new_file) THEN |
---|
[1492] | 327 | ! ++) Check number of month in file opened |
---|
| 328 | !************************************************************************************************** |
---|
| 329 | ierr = nf90_inq_dimid(ncid, 'TIME',dimid) |
---|
[1609] | 330 | CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps),"pb inq dim TIME" ) |
---|
[1492] | 331 | ! IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN |
---|
| 332 | IF (nbr_tsteps /= 12 ) THEN |
---|
[2317] | 333 | CALL abort_physic('get_aero_fromfile', 'not the right number of months in aerosol file read (should be 12 for the moment)' & |
---|
| 334 | ,1) |
---|
[1492] | 335 | ENDIF |
---|
[1150] | 336 | |
---|
[1179] | 337 | ! ++) Read the aerosol concentration month by month and concatenate to total variable varyear |
---|
| 338 | !**************************************************************************************** |
---|
| 339 | ! Get variable id |
---|
[1609] | 340 | CALL check_err( nf90_inq_varid(ncid, TRIM(varname), varid),"pb inq var "//TRIM(varname) ) |
---|
[1179] | 341 | |
---|
| 342 | ! Get the variable |
---|
[1609] | 343 | CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)),"pb get var "//TRIM(varname) ) |
---|
[1179] | 344 | |
---|
[1270] | 345 | ! ++) Read surface pression, 12 month in one variable |
---|
[1179] | 346 | !**************************************************************************************** |
---|
| 347 | ! Get variable id |
---|
[1609] | 348 | CALL check_err( nf90_inq_varid(ncid, "ps", varid),"pb inq var ps" ) |
---|
[1179] | 349 | ! Get the variable |
---|
[1609] | 350 | CALL check_err( nf90_get_var(ncid, varid, psurf_glo2D),"pb get var ps" ) |
---|
[1179] | 351 | |
---|
| 352 | ! ++) Read mass load, 12 month in one variable |
---|
| 353 | !**************************************************************************************** |
---|
| 354 | ! Get variable id |
---|
[1609] | 355 | CALL check_err( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) ,"pb inq var load_"//TRIM(varname)) |
---|
[1179] | 356 | ! Get the variable |
---|
[1609] | 357 | CALL check_err( nf90_get_var(ncid, varid, load_glo2D),"pb get var load_"//TRIM(varname) ) |
---|
[1179] | 358 | |
---|
| 359 | ! ++) Read ap |
---|
| 360 | !**************************************************************************************** |
---|
| 361 | ! Get variable id |
---|
[1609] | 362 | CALL check_err( nf90_inq_varid(ncid, "ap", varid),"pb inq var ap" ) |
---|
[1179] | 363 | ! Get the variable |
---|
[1609] | 364 | CALL check_err( nf90_get_var(ncid, varid, pt_ap),"pb get var ap" ) |
---|
[1150] | 365 | |
---|
[1179] | 366 | ! ++) Read b |
---|
| 367 | !**************************************************************************************** |
---|
| 368 | ! Get variable id |
---|
[1609] | 369 | CALL check_err( nf90_inq_varid(ncid, "b", varid),"pb inq var b" ) |
---|
[1179] | 370 | ! Get the variable |
---|
[1609] | 371 | CALL check_err( nf90_get_var(ncid, varid, pt_b),"pb get var b" ) |
---|
[1150] | 372 | |
---|
[1179] | 373 | ! ++) Read p0 : reference pressure |
---|
| 374 | !**************************************************************************************** |
---|
| 375 | ! Get variable id |
---|
[1609] | 376 | CALL check_err( nf90_inq_varid(ncid, "p0", varid),"pb inq var p0" ) |
---|
[1179] | 377 | ! Get the variable |
---|
[1609] | 378 | CALL check_err( nf90_get_var(ncid, varid, p0),"pb get var p0" ) |
---|
[1150] | 379 | |
---|
[524] | 380 | |
---|
[1179] | 381 | ELSE ! old file |
---|
[524] | 382 | |
---|
[1179] | 383 | ! ++) Read the aerosol concentration month by month and concatenate to total variable varyear |
---|
| 384 | !**************************************************************************************** |
---|
[1270] | 385 | DO imth=1, 12 |
---|
[1179] | 386 | IF (imth.EQ.1) THEN |
---|
| 387 | cvar=TRIM(varname)//'JAN' |
---|
| 388 | ELSE IF (imth.EQ.2) THEN |
---|
| 389 | cvar=TRIM(varname)//'FEB' |
---|
| 390 | ELSE IF (imth.EQ.3) THEN |
---|
| 391 | cvar=TRIM(varname)//'MAR' |
---|
| 392 | ELSE IF (imth.EQ.4) THEN |
---|
| 393 | cvar=TRIM(varname)//'APR' |
---|
| 394 | ELSE IF (imth.EQ.5) THEN |
---|
| 395 | cvar=TRIM(varname)//'MAY' |
---|
| 396 | ELSE IF (imth.EQ.6) THEN |
---|
| 397 | cvar=TRIM(varname)//'JUN' |
---|
| 398 | ELSE IF (imth.EQ.7) THEN |
---|
| 399 | cvar=TRIM(varname)//'JUL' |
---|
| 400 | ELSE IF (imth.EQ.8) THEN |
---|
| 401 | cvar=TRIM(varname)//'AUG' |
---|
| 402 | ELSE IF (imth.EQ.9) THEN |
---|
| 403 | cvar=TRIM(varname)//'SEP' |
---|
| 404 | ELSE IF (imth.EQ.10) THEN |
---|
| 405 | cvar=TRIM(varname)//'OCT' |
---|
| 406 | ELSE IF (imth.EQ.11) THEN |
---|
| 407 | cvar=TRIM(varname)//'NOV' |
---|
| 408 | ELSE IF (imth.EQ.12) THEN |
---|
| 409 | cvar=TRIM(varname)//'DEC' |
---|
| 410 | END IF |
---|
| 411 | |
---|
| 412 | ! Get variable id |
---|
[1609] | 413 | CALL check_err( nf90_inq_varid(ncid, TRIM(cvar), varid),"pb inq var "//TRIM(cvar) ) |
---|
[1179] | 414 | |
---|
| 415 | ! Get the variable |
---|
[1609] | 416 | CALL check_err( nf90_get_var(ncid, varid, varmth),"pb get var "//TRIM(cvar) ) |
---|
[1179] | 417 | |
---|
| 418 | ! Store in variable for the whole year |
---|
| 419 | varyear(:,:,:,imth)=varmth(:,:,:) |
---|
| 420 | |
---|
| 421 | END DO |
---|
[1150] | 422 | |
---|
[1179] | 423 | ! Putting dummy |
---|
| 424 | psurf_glo2D(:,:,:) = not_valid |
---|
| 425 | load_glo2D(:,:,:) = not_valid |
---|
| 426 | pt_ap(:) = not_valid |
---|
| 427 | pt_b(:) = not_valid |
---|
[524] | 428 | |
---|
[1179] | 429 | END IF |
---|
[524] | 430 | |
---|
[1179] | 431 | ! 4) Close file |
---|
| 432 | !**************************************************************************************** |
---|
[1609] | 433 | CALL check_err( nf90_close(ncid),"pb in close" ) |
---|
[1179] | 434 | |
---|
[524] | 435 | |
---|
[2346] | 436 | ! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo) |
---|
[1179] | 437 | !**************************************************************************************** |
---|
| 438 | ! Test if vertical levels have to be inversed |
---|
[782] | 439 | |
---|
[1179] | 440 | IF ((pt_b(1) < pt_b(klev_src)) .OR. .NOT. new_file) THEN |
---|
[1321] | 441 | ! WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' needs to be inverted' |
---|
| 442 | ! WRITE(lunout,*) 'before pt_ap = ', pt_ap |
---|
| 443 | ! WRITE(lunout,*) 'before pt_b = ', pt_b |
---|
[1179] | 444 | |
---|
| 445 | ! Inverse vertical levels for varyear |
---|
[1270] | 446 | DO imth=1, 12 |
---|
[1179] | 447 | varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly |
---|
| 448 | DO k=1, klev_src |
---|
[2346] | 449 | DO j=1, nbp_lat |
---|
| 450 | DO i=1, nbp_lon |
---|
[1179] | 451 | varyear(i,j,k,imth) = varmth(i,j,klev_src+1-k) |
---|
| 452 | END DO |
---|
| 453 | END DO |
---|
| 454 | END DO |
---|
| 455 | END DO |
---|
| 456 | |
---|
| 457 | ! Inverte vertical axes for pt_ap and pt_b |
---|
| 458 | varktmp(:) = pt_ap(:) |
---|
| 459 | DO k=1, klev_src |
---|
| 460 | pt_ap(k) = varktmp(klev_src+1-k) |
---|
| 461 | END DO |
---|
[782] | 462 | |
---|
[1179] | 463 | varktmp(:) = pt_b(:) |
---|
| 464 | DO k=1, klev_src |
---|
| 465 | pt_b(k) = varktmp(klev_src+1-k) |
---|
| 466 | END DO |
---|
| 467 | WRITE(lunout,*) 'after pt_ap = ', pt_ap |
---|
| 468 | WRITE(lunout,*) 'after pt_b = ', pt_b |
---|
[524] | 469 | |
---|
[1179] | 470 | ELSE |
---|
| 471 | WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done' |
---|
| 472 | WRITE(lunout,*) 'pt_ap = ', pt_ap |
---|
| 473 | WRITE(lunout,*) 'pt_b = ', pt_b |
---|
| 474 | END IF |
---|
[524] | 475 | |
---|
[1183] | 476 | ! - Invert latitudes if necessary |
---|
[1270] | 477 | DO imth=1, 12 |
---|
[1183] | 478 | IF (invert_lat) THEN |
---|
| 479 | |
---|
| 480 | ! Invert latitudes for the variable |
---|
| 481 | varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly |
---|
| 482 | DO k=1,klev_src |
---|
[2346] | 483 | DO j=1,nbp_lat |
---|
| 484 | DO i=1,nbp_lon |
---|
| 485 | varyear(i,j,k,imth) = varmth(i,nbp_lat+1-j,k) |
---|
[1183] | 486 | END DO |
---|
| 487 | END DO |
---|
| 488 | END DO |
---|
| 489 | |
---|
| 490 | ! Invert latitudes for surface pressure |
---|
| 491 | vartmp(:,:) = psurf_glo2D(:,:,imth) |
---|
[2346] | 492 | DO j=1,nbp_lat |
---|
| 493 | DO i=1,nbp_lon |
---|
| 494 | psurf_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j) |
---|
[1179] | 495 | END DO |
---|
| 496 | END DO |
---|
[1183] | 497 | |
---|
| 498 | ! Invert latitudes for the load |
---|
| 499 | vartmp(:,:) = load_glo2D(:,:,imth) |
---|
[2346] | 500 | DO j=1,nbp_lat |
---|
| 501 | DO i=1,nbp_lon |
---|
| 502 | load_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j) |
---|
[1183] | 503 | END DO |
---|
[1179] | 504 | END DO |
---|
[1183] | 505 | END IF ! invert_lat |
---|
| 506 | |
---|
[1179] | 507 | ! Do zonal mead at poles and distribut at whole first and last latitude |
---|
| 508 | DO k=1, klev_src |
---|
| 509 | npole=0. ! North pole, j=1 |
---|
[2346] | 510 | spole=0. ! South pole, j=nbp_lat |
---|
| 511 | DO i=1,nbp_lon |
---|
[1179] | 512 | npole = npole + varyear(i,1,k,imth) |
---|
[2346] | 513 | spole = spole + varyear(i,nbp_lat,k,imth) |
---|
[1179] | 514 | END DO |
---|
[2346] | 515 | npole = npole/REAL(nbp_lon) |
---|
| 516 | spole = spole/REAL(nbp_lon) |
---|
[1179] | 517 | varyear(:,1, k,imth) = npole |
---|
[2346] | 518 | varyear(:,nbp_lat,k,imth) = spole |
---|
[1179] | 519 | END DO |
---|
| 520 | END DO ! imth |
---|
[1183] | 521 | |
---|
[1270] | 522 | ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr) |
---|
[2311] | 523 | IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 3',1) |
---|
[1179] | 524 | |
---|
| 525 | ! Transform from 2D to 1D field |
---|
| 526 | CALL grid2Dto1D_glo(varyear,varyear_glo1D) |
---|
| 527 | CALL grid2Dto1D_glo(psurf_glo2D,psurf_glo1D) |
---|
| 528 | CALL grid2Dto1D_glo(load_glo2D,load_glo1D) |
---|
[1249] | 529 | |
---|
| 530 | ELSE |
---|
| 531 | ALLOCATE(varyear_glo1D(0,0,0)) |
---|
| 532 | END IF ! is_mpi_root .AND. is_omp_root |
---|
| 533 | |
---|
[1189] | 534 | !$OMP BARRIER |
---|
[1179] | 535 | |
---|
| 536 | ! 6) Distribute to all processes |
---|
| 537 | ! Scatter global field(klon_glo) to local process domain(klon) |
---|
| 538 | ! and distribute klev_src to all processes |
---|
| 539 | !**************************************************************************************** |
---|
[524] | 540 | |
---|
[1179] | 541 | ! Distribute klev_src |
---|
| 542 | CALL bcast(klev_src) |
---|
[524] | 543 | |
---|
[1179] | 544 | ! Allocate and distribute pt_ap and pt_b |
---|
| 545 | IF (.NOT. ASSOCIATED(pt_ap)) THEN ! if pt_ap is allocated also pt_b is allocated |
---|
| 546 | ALLOCATE(pt_ap(klev_src), pt_b(klev_src), stat=ierr) |
---|
[2311] | 547 | IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 4',1) |
---|
[1179] | 548 | END IF |
---|
| 549 | CALL bcast(pt_ap) |
---|
| 550 | CALL bcast(pt_b) |
---|
[524] | 551 | |
---|
[1179] | 552 | ! Allocate space for output pointer variable at local process |
---|
| 553 | IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year) |
---|
[1270] | 554 | ALLOCATE(pt_year(klon, klev_src, 12), stat=ierr) |
---|
[2311] | 555 | IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 5',1) |
---|
[524] | 556 | |
---|
[1179] | 557 | ! Scatter global field to local domain at local process |
---|
| 558 | CALL scatter(varyear_glo1D, pt_year) |
---|
[1270] | 559 | CALL scatter(psurf_glo1D, psurf_out) |
---|
| 560 | CALL scatter(load_glo1D, load_out) |
---|
[524] | 561 | |
---|
[1179] | 562 | ! 7) Test for negative values |
---|
| 563 | !**************************************************************************************** |
---|
| 564 | IF (MINVAL(pt_year) < 0.) THEN |
---|
| 565 | WRITE(lunout,*) 'Warning! Negative values read from file :', fname |
---|
| 566 | END IF |
---|
[524] | 567 | |
---|
[1179] | 568 | END SUBROUTINE get_aero_fromfile |
---|
[524] | 569 | |
---|
| 570 | |
---|
[1609] | 571 | SUBROUTINE check_err(status,text) |
---|
[1179] | 572 | USE netcdf |
---|
[2311] | 573 | USE print_control_mod, ONLY: lunout |
---|
[1179] | 574 | IMPLICIT NONE |
---|
[524] | 575 | |
---|
[1179] | 576 | INTEGER, INTENT (IN) :: status |
---|
[1609] | 577 | CHARACTER(len=*), INTENT (IN), OPTIONAL :: text |
---|
[524] | 578 | |
---|
[1179] | 579 | IF (status /= NF90_NOERR) THEN |
---|
[1609] | 580 | WRITE(lunout,*) 'Error in get_aero_fromfile, netcdf error code = ',status |
---|
| 581 | IF (PRESENT(text)) THEN |
---|
| 582 | WRITE(lunout,*) 'Error in get_aero_fromfile : ',text |
---|
| 583 | END IF |
---|
[2311] | 584 | CALL abort_physic('get_aero_fromfile',trim(nf90_strerror(status)),1) |
---|
[1179] | 585 | END IF |
---|
[524] | 586 | |
---|
[1179] | 587 | END SUBROUTINE check_err |
---|
| 588 | |
---|
| 589 | |
---|
| 590 | END MODULE readaerosol_mod |
---|