- Timestamp:
- Nov 20, 2009, 5:05:39 PM (15 years ago)
- Location:
- LMDZ4/branches/LMDZ4-dev/libf/phylmd
- Files:
-
- 2 edited
-
readaerosol.F90 (modified) (21 diffs)
-
readaerosol_interp.F90 (modified) (16 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol.F90
r1249 r1265 7 7 CONTAINS 8 8 9 SUBROUTINE readaerosol(name_aero, type, iyr_in, klev_src, pt_ap, pt_b, pt_out, p surf, load)9 SUBROUTINE readaerosol(name_aero, type, iyr_in, klev_src, pt_ap, pt_b, pt_out, pt_psurf, pt_load, nbr_tsteps) 10 10 11 11 !**************************************************************************************** … … 34 34 REAL, POINTER, DIMENSION(:) :: pt_ap ! Pointer for describing the vertical levels 35 35 REAL, POINTER, DIMENSION(:) :: pt_b ! Pointer for describing the vertical levels 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 36 REAL, POINTER, DIMENSION(:,:,:) :: pt_out ! The massvar distributions, DIMENSION(klon,klev_src, nbr_tsteps) 37 REAL, POINTER, DIMENSION(:,:) :: pt_psurf ! The massvar distributions, DIMENSION(klon, nbr_tsteps) 38 REAL, POINTER, DIMENSION(:,:) :: pt_load ! The massvar distributions, DIMENSION(klon, nbr_tsteps) 39 INTEGER :: nbr_tsteps ! number of timesteps in read file (12 or 14) 39 40 40 41 ! Local variables 41 CHARACTER(len=4) :: cyear42 REAL, POINTER, DIMENSION(:,:,:) :: pt_243 REAL, DIMENSION(klon,12):: psurf2, load244 REAL :: p0 ! Reference pressure45 INTEGER :: iyr1, iyr2, klev_src246 INTEGER :: it, k, i47 LOGICAL, PARAMETER :: lonlyone=.FALSE.42 CHARACTER(len=4) :: cyear 43 REAL, POINTER, DIMENSION(:,:,:) :: pt_2 44 REAL, POINTER, DIMENSION(:,:) :: psurf2, load2 45 REAL :: p0 ! Reference pressure 46 INTEGER :: iyr1, iyr2, klev_src2 47 INTEGER :: it, k, i 48 LOGICAL, PARAMETER :: lonlyone=.FALSE. 48 49 49 50 !**************************************************************************************** … … 56 57 !**************************************************************************************** 57 58 cyear='1980' 58 ! get_aero_fromfile returns pt_out allocated and initialized with data for 12month59 ! pt_out has dimensions (klon, klev_src, 12)60 CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, p surf, load)59 ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month 60 ! pt_out has dimensions (klon, klev_src, nbr_tsteps) 61 CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps) 61 62 62 63 … … 65 66 !**************************************************************************************** 66 67 cyear='.nat' 67 ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month 68 ! pt_out has dimensions (klon, klev_src, 12) 69 CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load) 68 ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month 69 ! pt_out has dimensions (klon, klev_src, nbr_tsteps) 70 CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps) 71 72 ELSE IF (type == 'annuel') THEN 73 ! Read and return data from scenario annual files 74 !**************************************************************************************** 75 WRITE(cyear,'(I4)') iyr_in 76 WRITE(lunout,*) 'get_aero 3 iyr_in=', iyr_in,' ',cyear 77 ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month 78 ! pt_out has dimensions (klon, klev_src, nbr_tsteps) 79 CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps) 80 70 81 71 82 ELSE IF (type == 'scenario') THEN … … 75 86 cyear='.nat' 76 87 WRITE(lunout,*) 'get_aero 1 iyr_in=', iyr_in,' ',cyear 77 ! get_aero_fromfile returns pt_out allocated and initialized with data for 12month78 ! pt_out has dimensions (klon, klev_src, 12)79 CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, p surf, load)88 ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month 89 ! pt_out has dimensions (klon, klev_src, nbr_tsteps) 90 CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps) 80 91 81 92 ELSE IF (iyr_in .GE. 2100) THEN 82 93 cyear='2100' 83 94 WRITE(lunout,*) 'get_aero 2 iyr_in=', iyr_in,' ',cyear 84 ! get_aero_fromfile returns pt_out allocated and initialized with data for 12month85 ! pt_out has dimensions (klon, klev_src, 12)86 CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, p surf, load)95 ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month 96 ! pt_out has dimensions (klon, klev_src, nbr_tsteps) 97 CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps) 87 98 88 99 ELSE … … 102 113 WRITE(cyear,'(I4)') iyr1 103 114 WRITE(lunout,*) 'get_aero 3 iyr_in=', iyr_in,' ',cyear 104 ! get_aero_fromfile returns pt_out allocated and initialized with data for 12month105 ! pt_out has dimensions (klon, klev_src, 12)106 CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, p surf, load)115 ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month 116 ! pt_out has dimensions (klon, klev_src, nbr_tsteps) 117 CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps) 107 118 108 119 ! If to read two decades: … … 114 125 115 126 NULLIFY(pt_2) 116 ! get_aero_fromfile returns pt_2 allocated and initialized with data for 12month117 ! pt_2 has dimensions (klon, klev_src, 12)118 CALL get_aero_fromfile(name_aero, cyear, klev_src2, pt_ap, pt_b, p0, pt_2, psurf2, load2 )127 ! get_aero_fromfile returns pt_2 allocated and initialized with data for nbr_tsteps month 128 ! pt_2 has dimensions (klon, klev_src, nbr_tsteps) 129 CALL get_aero_fromfile(name_aero, cyear, klev_src2, pt_ap, pt_b, p0, pt_2, psurf2, load2, nbr_tsteps) 119 130 ! Test for same number of vertical levels 120 131 IF (klev_src /= klev_src2) THEN … … 124 135 125 136 ! Linare interpolate to the actual year: 126 DO it=1, 12137 DO it=1,nbr_tsteps 127 138 DO k=1,klev_src 128 139 DO i = 1, klon … … 134 145 135 146 DO i = 1, klon 136 p surf(i,it) = &137 p surf(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &138 (p surf(i,it) - psurf2(i,it))139 140 load(i,it) = &141 load(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &142 ( load(i,it) - load2(i,it))147 pt_psurf(i,it) = & 148 pt_psurf(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * & 149 (pt_psurf(i,it) - psurf2(i,it)) 150 151 pt_load(i,it) = & 152 pt_load(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * & 153 (pt_load(i,it) - load2(i,it)) 143 154 END DO 144 155 END DO … … 146 157 ! Deallocate pt_2 no more needed 147 158 DEALLOCATE(pt_2) 159 DEALLOCATE(psurf2) 160 DEALLOCATE(load2) 148 161 149 162 END IF ! lonlyone … … 159 172 160 173 161 SUBROUTINE get_aero_fromfile(varname, cyr, klev_src, pt_ap, pt_b, p0, pt_year, p surf_out, load_out)162 !**************************************************************************************** 163 ! Read 12month aerosol from file and distribute to local process on physical grid.174 SUBROUTINE get_aero_fromfile(varname, cyr, klev_src, pt_ap, pt_b, p0, pt_year, pt_psurf_out, pt_load_out, nbr_tsteps) 175 !**************************************************************************************** 176 ! Read nbr_tsteps month aerosol from file and distribute to local process on physical grid. 164 177 ! Vertical levels, klev_src, may differ from model levels if new file format. 165 178 ! … … 197 210 REAL, POINTER, DIMENSION(:) :: pt_b ! Pointer for describing the vertical levels 198 211 REAL :: p0 ! Reference pressure value 199 REAL, POINTER, DIMENSION(:,:,:) :: pt_year ! Pointer-variabale from file, 12 month, grid : klon,klev_src 200 REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf_out ! Surface pression for 12 months 201 REAL, DIMENSION(klon,12), INTENT(OUT) :: load_out ! Aerosol mass load in each column 212 REAL, POINTER, DIMENSION(:,:,:) :: pt_year ! Pointer-variabale from file, nbr_tsteps month, grid : klon,klev_src 213 REAL, POINTER, DIMENSION(:,:) :: pt_psurf_out ! Surface pression for nbr_tsteps months 214 REAL, POINTER, DIMENSION(:,:) :: pt_load_out ! Aerosol mass load in each column 215 INTEGER :: nbr_tsteps ! number of month in file read 202 216 203 217 ! Local variables … … 209 223 REAL :: npole, spole 210 224 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: varmth 211 REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: varyear ! Global variable read from file, 12month212 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: varyear_glo1D !(klon_glo, klev_src, 12)225 REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: varyear ! Global variable read from file, nbr_tsteps month 226 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: varyear_glo1D !(klon_glo, klev_src, nbr_tsteps) 213 227 REAL, ALLOCATABLE, DIMENSION(:) :: varktmp 214 228 215 REAL, DIMENSION(iim,jjm+1,12) :: psurf_glo2D ! Surface pression for 12 months on dynamics global grid216 REAL, DIMENSION(klon_glo,12) :: psurf_glo1D ! -"- on physical global grid217 REAL, DIMENSION(iim,jjm+1,12) :: load_glo2D ! Load for 12 months on dynamics global grid218 REAL, DIMENSION(klon_glo,12) :: load_glo1D ! -"- on physical global grid229 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: psurf_glo2D 230 REAL, ALLOCATABLE, DIMENSION(:,:) :: psurf_glo1D 231 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: load_glo2D 232 REAL, ALLOCATABLE, DIMENSION(:,:) :: load_glo1D 219 233 REAL, DIMENSION(iim,jjm+1) :: vartmp 220 234 REAL, DIMENSION(iim) :: lon_src ! longitudes in file … … 273 287 CALL abort_gcm('get_aero_fromfile', 'latitudes do not correspond between file and model',1) 274 288 END IF 289 290 ! 1.5) Check number of month in file opened 291 ! 292 !************************************************************************************************** 293 ierr = nf90_inq_dimid(ncid, 'TIME',dimid) 294 ! ierr = nf90_inq_dimlen(ncid, dimid, nbr_tsteps) 295 CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps) ) 296 IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN 297 CALL abort_gcm('get_aero_fromfile', 'not the right number of months in aerosol file read',1) 298 ENDIF 299 write(lunout,*)'get_aero_fromfile: ', nbr_tsteps,' months to read from aerosols file' 300 call bcast(nbr_tsteps) 301 302 ! 1.6) Allocation of some variables once the number of months in the file has been determined 303 ! 304 !************************************************************************************************** 305 ALLOCATE(psurf_glo2D(iim, jjm+1, nbr_tsteps), load_glo2D(iim, jjm+1, nbr_tsteps), stat=ierr) 306 IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 0.5',1) 307 ALLOCATE(psurf_glo1D(klon_glo, nbr_tsteps), load_glo1D(klon_glo, nbr_tsteps), stat=ierr) 308 IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 0.5',1) 309 275 310 276 311 ! 2) Check if old or new file is avalabale. … … 301 336 302 337 ! Allocate variables depending on the number of vertical levels 303 ALLOCATE(varmth(iim, jjm+1, klev_src), varyear(iim, jjm+1, klev_src, 12), stat=ierr)338 ALLOCATE(varmth(iim, jjm+1, klev_src), varyear(iim, jjm+1, klev_src, nbr_tsteps), stat=ierr) 304 339 IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 1',1) 305 340 … … 323 358 CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)) ) 324 359 325 ! ++) Read surface pression, 12month in one variable360 ! ++) Read surface pression, nbr_tsteps month in one variable 326 361 !**************************************************************************************** 327 362 ! Get variable id … … 363 398 ! ++) Read the aerosol concentration month by month and concatenate to total variable varyear 364 399 !**************************************************************************************** 365 DO imth=1, 12400 DO imth=1, nbr_tsteps 366 401 IF (imth.EQ.1) THEN 367 402 cvar=TRIM(varname)//'JAN' … … 424 459 425 460 ! Inverse vertical levels for varyear 426 DO imth=1, 12461 DO imth=1, nbr_tsteps 427 462 varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly 428 463 DO k=1, klev_src … … 455 490 456 491 ! - Invert latitudes if necessary 457 DO imth=1, 12492 DO imth=1, nbr_tsteps 458 493 IF (invert_lat) THEN 459 494 … … 500 535 END DO ! imth 501 536 502 ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr)537 ALLOCATE(varyear_glo1D(klon_glo, klev_src, nbr_tsteps), stat=ierr) 503 538 IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 3',1) 504 539 … … 532 567 ! Allocate space for output pointer variable at local process 533 568 IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year) 534 ALLOCATE(pt_year(klon, klev_src, 12), stat=ierr) 535 IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 5',1) 569 write(lunout,*)'get_aero_fromfile: klon, ...',klon, klev_src, nbr_tsteps 570 ALLOCATE(pt_year(klon, klev_src, nbr_tsteps), stat=ierr) 571 ! IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 5.1',1) 572 IF (ASSOCIATED(pt_psurf_out)) DEALLOCATE(pt_psurf_out) 573 ALLOCATE(pt_psurf_out(klon, nbr_tsteps), stat=ierr) 574 ! IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 5.2',1) 575 IF (ASSOCIATED(pt_load_out)) DEALLOCATE(pt_load_out) 576 ALLOCATE(pt_load_out(klon, nbr_tsteps), stat=ierr) 577 ! IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 5.3',1) 536 578 537 579 ! Scatter global field to local domain at local process 538 580 CALL scatter(varyear_glo1D, pt_year) 539 CALL scatter(psurf_glo1D, p surf_out)540 CALL scatter(load_glo1D, load_out)581 CALL scatter(psurf_glo1D, pt_psurf_out) 582 CALL scatter(load_glo1D, pt_load_out) 541 583 542 584 ! 7) Test for negative values -
LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol_interp.F90
r1249 r1265 13 13 ! 4) Test for negative mass values 14 14 15 USE ioipsl 15 16 USE dimphy, ONLY : klev,klon 16 17 USE mod_phys_lmdz_para, ONLY : mpi_rank … … 50 51 INTEGER :: i, k, ierr 51 52 INTEGER :: iday, iyr, lmt_pas 52 INTEGER :: im, day1, day2, im2 53 ! INTEGER :: im, day1, day2, im2 54 INTEGER :: im, im2 55 REAL :: day1, day2 53 56 INTEGER :: pi_klev_src ! Only for testing purpose 54 57 INTEGER, SAVE :: klev_src ! Number of vertical levles in source field … … 76 79 77 80 REAL, DIMENSION(:,:,:), POINTER :: pt_tmp ! Pointer allocated in readaerosol 81 REAL, DIMENSION(:,:), POINTER :: pt_tmp_surf ! Pointer allocated in readaerosol 82 REAL, DIMENSION(:,:), POINTER :: pt_tmp_load ! Pointer allocated in readaerosol 83 REAL, DIMENSION(:,:), POINTER :: pt_tmp_pi_surf ! Pointer allocated in readaerosol 84 REAL, DIMENSION(:,:), POINTER :: pt_tmp_pi_load ! Pointer allocated in readaerosol 78 85 REAL, POINTER, DIMENSION(:), SAVE :: pt_ap, pt_b ! Pointer for describing the vertical levels 79 86 !$OMP THREADPRIVATE(pt_ap, pt_b) 80 87 INTEGER, SAVE :: nbr_tsteps ! number of time steps in file read 88 REAL, DIMENSION(14), SAVE :: month_len, month_start, month_mid 89 !$OMP THREADPRIVATE(month_len, month_start, month_mid) 90 REAL :: jDay 91 92 !$OMP THREADPRIVATE(nbr_tsteps) 81 93 LOGICAL :: lnewday ! Indicates if first time step at a new day 82 94 LOGICAL :: OLDNEWDAY … … 98 110 99 111 ! Use phys_cal_mod 100 !iday= day_cur101 !iyr = year_cur102 !im = mth_cur112 ! iday= day_cur 113 ! iyr = year_cur 114 ! im = mth_cur 103 115 104 116 iday = INT(r_day) … … 107 119 iyr = iyr + annee_ref ! year of the run 108 120 im = iday/30 +1 ! the actual month 121 CALL ymds2ju(iyr, im, iday, 0., jDay) 109 122 110 123 IF(MOD(itap-1,NINT(86400./pdtphys)) == 0)THEN … … 140 153 IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 2',1) 141 154 142 ALLOCATE( psurf_year(klon, 12, naero_spc), pi_psurf_year(klon, 12, naero_spc), stat=ierr)143 IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 3',1)144 145 ALLOCATE( load_year(klon, 12, naero_spc), pi_load_year(klon, 12, naero_spc), stat=ierr)146 IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 4',1)155 ! ALLOCATE( psurf_year(klon, 12, naero_spc), pi_psurf_year(klon, 12, naero_spc), stat=ierr) 156 ! IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 3',1) 157 ! 158 ! ALLOCATE( load_year(klon, 12, naero_spc), pi_load_year(klon, 12, naero_spc), stat=ierr) 159 ! IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 4',1) 147 160 148 161 lnewday=.TRUE. … … 159 172 IF ( (first .OR. iday==0) .AND. lnewday ) THEN 160 173 NULLIFY(pt_tmp) 174 NULLIFY(pt_tmp_surf) 175 NULLIFY(pt_tmp_pi_surf) 176 NULLIFY(pt_tmp_load) 177 NULLIFY(pt_tmp_pi_load) 161 178 162 179 ! Reading values corresponding to the closest year taking into count the choice of aer_type. 163 180 ! For aer_type=scenario interpolation between 2 data sets is done in readaerosol. 164 181 CALL readaerosol(name_aero(id_aero), aer_type, iyr, klev_src, pt_ap, pt_b, pt_tmp, & 165 p surf_year(:,:,id_aero), load_year(:,:,id_aero))182 pt_tmp_surf,pt_tmp_load,nbr_tsteps) 166 183 IF (.NOT. ALLOCATED(var_year)) THEN 167 ALLOCATE(var_year(klon, klev_src, 12, naero_spc), stat=ierr) 168 IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 5',1) 184 ALLOCATE(var_year(klon, klev_src, nbr_tsteps, naero_spc), stat=ierr) 185 IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 5.1',1) 186 ALLOCATE( psurf_year(klon, nbr_tsteps, naero_spc), load_year(klon, nbr_tsteps, naero_spc), stat=ierr) 187 IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 5.2',1) 169 188 END IF 170 189 var_year(:,:,:,id_aero) = pt_tmp(:,:,:) 190 psurf_year(:,:,id_aero) = pt_tmp_surf(:,:) 191 load_year(:,:,id_aero) = pt_tmp_load(:,:) 171 192 172 193 ! Reading values corresponding to the preindustrial concentrations. 173 194 CALL readaerosol(name_aero(id_aero), 'preind', iyr, pi_klev_src, pt_ap, pt_b, pt_tmp, & 174 p i_psurf_year(:,:,id_aero), pi_load_year(:,:,id_aero))195 pt_tmp_surf,pt_tmp_load,nbr_tsteps) 175 196 176 197 ! klev_src must be the same in both files. … … 183 204 184 205 IF (.NOT. ALLOCATED(pi_var_year)) THEN 185 ALLOCATE(pi_var_year(klon, klev_src, 12, naero_spc), stat=ierr) 186 IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 6',1) 206 ALLOCATE(pi_var_year(klon, klev_src, nbr_tsteps, naero_spc), stat=ierr) 207 IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 6.1',1) 208 ALLOCATE( pi_psurf_year(klon, nbr_tsteps, naero_spc), pi_load_year(klon, nbr_tsteps, naero_spc), stat=ierr) 209 IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 6.2',1) 187 210 END IF 188 211 pi_var_year(:,:,:,id_aero) = pt_tmp(:,:,:) 212 pi_psurf_year(:,:,id_aero) = pt_tmp_surf(:,:) 213 pi_load_year(:,:,id_aero) = pt_tmp_load(:,:) 189 214 190 215 IF (debug) THEN 191 CALL writefield_phy('var_year_ jan',var_year(:,:,1,id_aero),klev_src)192 CALL writefield_phy('var_year_ dec',var_year(:,:,12,id_aero),klev_src)216 CALL writefield_phy('var_year_first',var_year(:,:,1,id_aero),klev_src) 217 CALL writefield_phy('var_year_last',var_year(:,:,nbr_tsteps,id_aero),klev_src) 193 218 CALL writefield_phy('psurf_src',psurf_year(:,:,id_aero),1) 194 219 CALL writefield_phy('pi_psurf_src',pi_psurf_year(:,:,id_aero),1) … … 199 224 ! Pointer no more useful, deallocate. 200 225 DEALLOCATE(pt_tmp) 226 DEALLOCATE(pt_tmp_surf) 227 DEALLOCATE(pt_tmp_load) 201 228 202 229 ! Test if vertical interpolation will be needed. … … 220 247 END IF 221 248 249 ! Calendar initialisation 250 ! 251 DO i = 2, 13 252 month_len(i) = float(ioget_mon_len(year_cur, i-1)) 253 CALL ymds2ju(year_cur, i-1, 1, 0.0, month_start(i)) 254 ENDDO 255 month_len(1) = float(ioget_mon_len(year_cur-1, 12)) 256 CALL ymds2ju(year_cur-1, 12, 1, 0.0, month_start(1)) 257 month_len(14) = float(ioget_mon_len(year_cur+1, 1)) 258 CALL ymds2ju(year_cur+1, 1, 1, 0.0, month_start(14)) 259 month_mid(:) = month_start (:) + month_len(:)/2. 260 261 write(55,*)'month_len = ',month_len 262 write(55,*)'month_start = ',month_start 263 write(55,*)'month_mid = ',month_mid 264 222 265 END IF ! IF ( (first .OR. iday==0) .AND. lnewday ) THEN 223 266 … … 233 276 ! 234 277 !**************************************************************************************** 235 ! Find which months and days to use for time interpolation 236 IF (iday < im*30-15) THEN 237 ! in the first half of the month use month before and actual month 238 im2=im-1 239 day2 = im2*30-15 240 day1 = im2*30+15 241 IF (im2 <= 0) THEN 242 ! the month is january, thus the month before december 243 im2=12 244 END IF 278 ! Find which months and days to use for time interpolation 279 IF (nbr_tsteps == 12) then 280 IF (jDay < month_mid(im+1)) THEN 281 im2=im-1 282 day2 = month_mid(im2+1) 283 day1 = month_mid(im+1) 284 IF (im2 <= 0) THEN 285 ! the month is january, thus the month before december 286 im2=12 287 END IF 288 ELSE 289 ! the second half of the month 290 im2=im+1 291 day2 = month_mid(im+1) 292 day1 = month_mid(im2+1) 293 IF (im2 > 12) THEN 294 ! the month is december, the following thus january 295 im2=1 296 ENDIF 297 END IF 298 ELSE IF (nbr_tsteps == 14) then 299 im = im + 1 300 IF (jDay < month_mid(im)) THEN 301 ! in the first half of the month use month before and actual month 302 im2=im-1 303 day2 = month_mid(im2) 304 day1 = month_mid(im) 305 ELSE 306 ! the second half of the month 307 im2=im+1 308 day2 = month_mid(im) 309 day1 = month_mid(im2) 310 END IF 245 311 ELSE 246 ! the second half of the month 247 im2=im+1 248 IF (im2 > 12) THEN 249 ! the month is december, the following thus january 250 im2=1 251 ENDIF 252 day2 = im*30-15 253 day1 = im*30+15 254 END IF 255 312 CALL abort_gcm('readaerosol_interp', 'number of months undefined',1) 313 ENDIF 314 315 ! ! Find which months and days to use for time interpolation 316 ! IF (iday < im*30-15) THEN 317 ! ! in the first half of the month use month before and actual month 318 ! im2=im-1 319 ! day2 = im2*30-15 320 ! day1 = im2*30+15 321 ! IF (im2 <= 0) THEN 322 ! ! the month is january, thus the month before december 323 ! im2=12 324 ! END IF 325 ! ELSE 326 ! ! the second half of the month 327 ! im2=im+1 328 ! IF (im2 > 12) THEN 329 ! ! the month is december, the following thus january 330 ! im2=1 331 ! ENDIF 332 ! day2 = im*30-15 333 ! day1 = im*30+15 334 ! END IF 335 ! jDay = jDay+1 336 ! write(55,*)'iday, jDay, im, im2, day1, day2 = ',iday, jDay, im, im2, day1, day2, nbr_tsteps 337 256 338 ! Time interpolation, still on vertical source grid 257 339 ALLOCATE(tmp1(klon,klev_src), tmp2(klon,klev_src),stat=ierr) … … 265 347 DO i=1,klon 266 348 tmp1(i,k) = & 267 var_year(i,k,im2,id_aero) - FLOAT( iday-day2)/FLOAT(day1-day2) * &349 var_year(i,k,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * & 268 350 (var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero)) 269 351 270 352 tmp2(i,k) = & 271 pi_var_year(i,k,im2,id_aero) - FLOAT( iday-day2)/FLOAT(day1-day2) * &353 pi_var_year(i,k,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * & 272 354 (pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero)) 273 355 END DO … … 277 359 DO i=1,klon 278 360 psurf_day(i) = & 279 psurf_year(i,im2,id_aero) - FLOAT( iday-day2)/FLOAT(day1-day2) * &361 psurf_year(i,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * & 280 362 (psurf_year(i,im2,id_aero) - psurf_year(i,im,id_aero)) 281 363 282 364 pi_psurf_day(i) = & 283 pi_psurf_year(i,im2,id_aero) - FLOAT( iday-day2)/FLOAT(day1-day2) * &365 pi_psurf_year(i,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * & 284 366 (pi_psurf_year(i,im2,id_aero) - pi_psurf_year(i,im,id_aero)) 285 367 END DO … … 288 370 DO i=1,klon 289 371 load_src(i) = & 290 load_year(i,im2,id_aero) - FLOAT( iday-day2)/FLOAT(day1-day2) * &372 load_year(i,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * & 291 373 (load_year(i,im2,id_aero) - load_year(i,im,id_aero)) 292 374 293 375 pi_load_src(i) = & 294 pi_load_year(i,im2,id_aero) - FLOAT( iday-day2)/FLOAT(day1-day2) * &376 pi_load_year(i,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * & 295 377 (pi_load_year(i,im2,id_aero) - pi_load_year(i,im,id_aero)) 296 378 END DO … … 447 529 ! Test for var_day 448 530 IF (var_day(i,k,id_aero) < 0.) THEN 449 IF ( iday-day2 < 0.) WRITE(lunout,*) 'iday-day2=',iday-day2531 IF (jDay-day2 < 0.) WRITE(lunout,*) 'jDay-day2=',jDay-day2 450 532 IF (var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero) < 0.) THEN 451 533 WRITE(lunout,*) trim(name_aero(id_aero)),'(i,k,im2)-', & … … 466 548 ! Test for pi_var_day 467 549 IF (pi_var_day(i,k,id_aero) < 0.) THEN 468 IF ( iday-day2 < 0.) WRITE(lunout,*) 'iday-day2=',iday-day2550 IF (jDay-day2 < 0.) WRITE(lunout,*) 'jDay-day2=',jDay-day2 469 551 IF (pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero) < 0.) THEN 470 552 WRITE(lunout,*) trim(name_aero(id_aero)),'(i,k,im2)-', &
Note: See TracChangeset
for help on using the changeset viewer.
