- Timestamp:
- Nov 20, 2009, 5:05:39 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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.
