| Line |
Branch |
Exec |
Source |
| 1 |
|
|
! $Id$ |
| 2 |
|
|
! |
| 3 |
|
✗ |
SUBROUTINE readaerosol_interp(id_aero, itap, pdtphys, r_day, first, pplay, paprs, t_seri, mass_out, pi_mass_out, load_src) |
| 4 |
|
|
! |
| 5 |
|
|
! This routine will return the mass concentration at actual day(mass_out) and |
| 6 |
|
|
! the pre-industrial values(pi_mass_out) for aerosol corresponding to "id_aero". |
| 7 |
|
|
! The mass concentrations for all aerosols are saved in this routine but each |
| 8 |
|
|
! call to this routine only treats the aerosol "id_aero". |
| 9 |
|
|
! |
| 10 |
|
|
! 1) Read in data for the whole year, only at first time step |
| 11 |
|
|
! 2) Interpolate to the actual day, only at new day |
| 12 |
|
|
! 3) Interpolate to the model vertical grid (target grid), only at new day |
| 13 |
|
|
! 4) Test for negative mass values |
| 14 |
|
|
|
| 15 |
|
|
USE ioipsl |
| 16 |
|
|
USE dimphy, ONLY : klev,klon |
| 17 |
|
|
USE mod_phys_lmdz_para, ONLY : mpi_rank |
| 18 |
|
|
USE readaerosol_mod |
| 19 |
|
|
USE aero_mod, ONLY : naero_spc, name_aero |
| 20 |
|
|
USE write_field_phy |
| 21 |
|
|
USE phys_cal_mod |
| 22 |
|
|
USE pres2lev_mod |
| 23 |
|
|
USE print_control_mod, ONLY: lunout |
| 24 |
|
|
|
| 25 |
|
|
IMPLICIT NONE |
| 26 |
|
|
|
| 27 |
|
|
INCLUDE "YOMCST.h" |
| 28 |
|
|
INCLUDE "chem.h" |
| 29 |
|
|
INCLUDE "clesphys.h" |
| 30 |
|
|
|
| 31 |
|
|
! |
| 32 |
|
|
! Input: |
| 33 |
|
|
!**************************************************************************************** |
| 34 |
|
|
INTEGER, INTENT(IN) :: id_aero! Identity number for the aerosol to treat |
| 35 |
|
|
INTEGER, INTENT(IN) :: itap ! Physic step count |
| 36 |
|
|
REAL, INTENT(IN) :: pdtphys! Physic day step |
| 37 |
|
|
REAL, INTENT(IN) :: r_day ! Day of integration |
| 38 |
|
|
LOGICAL, INTENT(IN) :: first ! First model timestep |
| 39 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! pression at model mid-layers |
| 40 |
|
|
REAL, DIMENSION(klon,klev+1),INTENT(IN):: paprs ! pression between model layers |
| 41 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: t_seri ! air temperature |
| 42 |
|
|
! |
| 43 |
|
|
! Output: |
| 44 |
|
|
!**************************************************************************************** |
| 45 |
|
|
REAL, INTENT(OUT) :: mass_out(klon,klev) ! Mass of aerosol (monthly mean data,from file) [ug AIBCM/m3] |
| 46 |
|
|
REAL, INTENT(OUT) :: pi_mass_out(klon,klev) ! Mass of preindustrial aerosol (monthly mean data,from file) [ug AIBCM/m3] |
| 47 |
|
|
REAL, INTENT(OUT) :: load_src(klon) ! Load of aerosol (monthly mean data,from file) [kg/m3] |
| 48 |
|
|
! |
| 49 |
|
|
! Local Variables: |
| 50 |
|
|
!**************************************************************************************** |
| 51 |
|
|
INTEGER :: i, k, ierr |
| 52 |
|
|
INTEGER :: iday, iyr, lmt_pas |
| 53 |
|
|
! INTEGER :: im, day1, day2, im2 |
| 54 |
|
|
INTEGER :: im, im2 |
| 55 |
|
|
REAL :: day1, day2 |
| 56 |
|
|
INTEGER :: pi_klev_src ! Only for testing purpose |
| 57 |
|
|
INTEGER, SAVE :: klev_src ! Number of vertical levles in source field |
| 58 |
|
|
!$OMP THREADPRIVATE(klev_src) |
| 59 |
|
|
|
| 60 |
|
|
REAL :: zrho ! Air density [kg/m3] |
| 61 |
|
|
REAL :: volm ! Volyme de melange [kg/kg] |
| 62 |
|
✗ |
REAL, DIMENSION(klon) :: psurf_day, pi_psurf_day |
| 63 |
|
✗ |
REAL, DIMENSION(klon) :: pi_load_src ! Mass load at source grid |
| 64 |
|
✗ |
REAL, DIMENSION(klon) :: load_tgt, load_tgt_test |
| 65 |
|
✗ |
REAL, DIMENSION(klon,klev) :: delp ! pressure difference in each model layer |
| 66 |
|
|
|
| 67 |
|
|
REAL, ALLOCATABLE, DIMENSION(:,:) :: pplay_src ! pression mid-layer at source levels |
| 68 |
|
|
REAL, ALLOCATABLE, DIMENSION(:,:) :: tmp1, tmp2 ! Temporary variables |
| 69 |
|
|
REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: var_year ! VAR in right dimension for the total year |
| 70 |
|
|
REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: pi_var_year ! pre-industrial VAR, -"- |
| 71 |
|
|
!$OMP THREADPRIVATE(var_year,pi_var_year) |
| 72 |
|
|
REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: var_day ! VAR interpolated to the actual day and model grid |
| 73 |
|
|
REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: pi_var_day ! pre-industrial VAR, -"- |
| 74 |
|
|
!$OMP THREADPRIVATE(var_day,pi_var_day) |
| 75 |
|
|
REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: psurf_year, pi_psurf_year ! surface pressure for the total year |
| 76 |
|
|
!$OMP THREADPRIVATE(psurf_year, pi_psurf_year) |
| 77 |
|
|
REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: load_year, pi_load_year ! load in the column for the total year |
| 78 |
|
|
!$OMP THREADPRIVATE(load_year, pi_load_year) |
| 79 |
|
|
|
| 80 |
|
|
REAL, DIMENSION(:,:,:), POINTER :: pt_tmp ! Pointer allocated in readaerosol |
| 81 |
|
|
REAL, POINTER, DIMENSION(:), SAVE :: pt_ap, pt_b ! Pointer for describing the vertical levels |
| 82 |
|
|
!$OMP THREADPRIVATE(pt_ap, pt_b) |
| 83 |
|
|
INTEGER, SAVE :: nbr_tsteps ! number of time steps in file read |
| 84 |
|
|
REAL, DIMENSION(14), SAVE :: month_len, month_start, month_mid |
| 85 |
|
|
!$OMP THREADPRIVATE(nbr_tsteps, month_len, month_start, month_mid) |
| 86 |
|
|
REAL :: jDay |
| 87 |
|
|
|
| 88 |
|
|
LOGICAL :: lnewday ! Indicates if first time step at a new day |
| 89 |
|
|
LOGICAL :: OLDNEWDAY |
| 90 |
|
|
LOGICAL,SAVE :: vert_interp ! Indicates if vertical interpolation will be done |
| 91 |
|
|
LOGICAL,SAVE :: debug=.FALSE.! Debugging in this subroutine |
| 92 |
|
|
!$OMP THREADPRIVATE(vert_interp, debug) |
| 93 |
|
|
CHARACTER(len=8) :: type |
| 94 |
|
|
CHARACTER(len=8) :: filename |
| 95 |
|
|
|
| 96 |
|
|
|
| 97 |
|
|
!**************************************************************************************** |
| 98 |
|
|
! Initialization |
| 99 |
|
|
! |
| 100 |
|
|
!**************************************************************************************** |
| 101 |
|
|
|
| 102 |
|
|
! Calculation to find if it is a new day |
| 103 |
|
|
|
| 104 |
|
✗ |
IF(mpi_rank == 0 .AND. debug )then |
| 105 |
|
✗ |
PRINT*,'CONTROL PANEL REGARDING TIME STEPING' |
| 106 |
|
|
ENDIF |
| 107 |
|
|
|
| 108 |
|
|
! Use phys_cal_mod |
| 109 |
|
✗ |
iday= day_cur |
| 110 |
|
✗ |
iyr = year_cur |
| 111 |
|
✗ |
im = mth_cur |
| 112 |
|
|
|
| 113 |
|
|
! iday = INT(r_day) |
| 114 |
|
|
! iyr = iday/360 |
| 115 |
|
|
! iday = iday-iyr*360 ! day of the actual year |
| 116 |
|
|
! iyr = iyr + annee_ref ! year of the run |
| 117 |
|
|
! im = iday/30 +1 ! the actual month |
| 118 |
|
✗ |
CALL ymds2ju(iyr, im, iday, 0., jDay) |
| 119 |
|
|
! CALL ymds2ju(iyr, im, iday-(im-1)*30, 0., jDay) |
| 120 |
|
|
|
| 121 |
|
|
|
| 122 |
|
✗ |
IF(MOD(itap-1,NINT(86400./pdtphys)) == 0)THEN |
| 123 |
|
✗ |
lnewday=.TRUE. |
| 124 |
|
|
ELSE |
| 125 |
|
✗ |
lnewday=.FALSE. |
| 126 |
|
|
ENDIF |
| 127 |
|
|
|
| 128 |
|
✗ |
IF(mpi_rank == 0 .AND. debug)then |
| 129 |
|
|
! 0.02 is about 0.5/24, namly less than half an hour |
| 130 |
|
✗ |
OLDNEWDAY = (r_day-REAL(iday) < 0.02) |
| 131 |
|
|
! Once per day, update aerosol fields |
| 132 |
|
✗ |
lmt_pas = NINT(86400./pdtphys) |
| 133 |
|
✗ |
PRINT*,'r_day-REAL(iday) =',r_day-REAL(iday) |
| 134 |
|
✗ |
PRINT*,'itap =',itap |
| 135 |
|
✗ |
PRINT*,'pdtphys =',pdtphys |
| 136 |
|
✗ |
PRINT*,'lmt_pas =',lmt_pas |
| 137 |
|
✗ |
PRINT*,'iday =',iday |
| 138 |
|
✗ |
PRINT*,'r_day =',r_day |
| 139 |
|
✗ |
PRINT*,'day_cur =',day_cur |
| 140 |
|
✗ |
PRINT*,'mth_cur =',mth_cur |
| 141 |
|
✗ |
PRINT*,'year_cur =',year_cur |
| 142 |
|
✗ |
PRINT*,'NINT(86400./pdtphys) =',NINT(86400./pdtphys) |
| 143 |
|
✗ |
PRINT*,'MOD(0,1) =',MOD(0,1) |
| 144 |
|
✗ |
PRINT*,'lnewday =',lnewday |
| 145 |
|
✗ |
PRINT*,'OLDNEWDAY =',OLDNEWDAY |
| 146 |
|
|
ENDIF |
| 147 |
|
|
|
| 148 |
|
✗ |
IF (.NOT. ALLOCATED(var_day)) THEN |
| 149 |
|
✗ |
ALLOCATE( var_day(klon, klev, naero_spc), stat=ierr) |
| 150 |
|
✗ |
IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 1',1) |
| 151 |
|
✗ |
ALLOCATE( pi_var_day(klon, klev, naero_spc), stat=ierr) |
| 152 |
|
✗ |
IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 2',1) |
| 153 |
|
|
|
| 154 |
|
✗ |
ALLOCATE( psurf_year(klon, 12, naero_spc), pi_psurf_year(klon, 12, naero_spc), stat=ierr) |
| 155 |
|
✗ |
IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 3',1) |
| 156 |
|
|
|
| 157 |
|
✗ |
ALLOCATE( load_year(klon, 12, naero_spc), pi_load_year(klon, 12, naero_spc), stat=ierr) |
| 158 |
|
✗ |
IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 4',1) |
| 159 |
|
|
|
| 160 |
|
✗ |
lnewday=.TRUE. |
| 161 |
|
|
|
| 162 |
|
✗ |
NULLIFY(pt_ap) |
| 163 |
|
✗ |
NULLIFY(pt_b) |
| 164 |
|
|
ENDIF |
| 165 |
|
|
|
| 166 |
|
|
!**************************************************************************************** |
| 167 |
|
|
! 1) Read in data : corresponding to the actual year and preindustrial data. |
| 168 |
|
|
! Only for the first day of the year. |
| 169 |
|
|
! |
| 170 |
|
|
!**************************************************************************************** |
| 171 |
|
✗ |
IF ( (first .OR. iday==0) .AND. lnewday ) THEN |
| 172 |
|
✗ |
NULLIFY(pt_tmp) |
| 173 |
|
|
|
| 174 |
|
|
! Reading values corresponding to the closest year taking into count the choice of aer_type. |
| 175 |
|
|
! For aer_type=scenario interpolation between 2 data sets is done in readaerosol. |
| 176 |
|
|
! If aer_type=mix1, mix2 or mix3, the run type and file name depends on the aerosol. |
| 177 |
|
✗ |
IF (aer_type=='preind' .OR. aer_type=='actuel' .OR. aer_type=='annuel' .OR. aer_type=='scenario') THEN |
| 178 |
|
|
! Standard case |
| 179 |
|
✗ |
filename='aerosols' |
| 180 |
|
✗ |
type=aer_type |
| 181 |
|
✗ |
ELSE IF (aer_type == 'mix1') THEN |
| 182 |
|
|
! Special case using a mix of decenal sulfate file and annual aerosols(all aerosols except sulfate) |
| 183 |
|
✗ |
IF (name_aero(id_aero) == 'SO4') THEN |
| 184 |
|
✗ |
filename='so4.run ' |
| 185 |
|
✗ |
type='scenario' |
| 186 |
|
|
ELSE |
| 187 |
|
✗ |
filename='aerosols' |
| 188 |
|
✗ |
type='annuel' |
| 189 |
|
|
ENDIF |
| 190 |
|
✗ |
ELSE IF (aer_type == 'mix2') THEN |
| 191 |
|
|
! Special case using a mix of decenal sulfate file and natrual aerosols |
| 192 |
|
✗ |
IF (name_aero(id_aero) == 'SO4') THEN |
| 193 |
|
✗ |
filename='so4.run ' |
| 194 |
|
✗ |
type='scenario' |
| 195 |
|
|
ELSE |
| 196 |
|
✗ |
filename='aerosols' |
| 197 |
|
✗ |
type='preind' |
| 198 |
|
|
ENDIF |
| 199 |
|
✗ |
ELSE IF (aer_type == 'mix3') THEN |
| 200 |
|
|
! Special case using a mix of annual sulfate file and natrual aerosols |
| 201 |
|
✗ |
IF (name_aero(id_aero) == 'SO4') THEN |
| 202 |
|
✗ |
filename='aerosols' |
| 203 |
|
✗ |
type='annuel' |
| 204 |
|
|
ELSE |
| 205 |
|
✗ |
filename='aerosols' |
| 206 |
|
✗ |
type='preind' |
| 207 |
|
|
ENDIF |
| 208 |
|
|
ELSE |
| 209 |
|
✗ |
CALL abort_physic('readaerosol_interp', 'this aer_type not supported',1) |
| 210 |
|
|
ENDIF |
| 211 |
|
|
|
| 212 |
|
|
CALL readaerosol(name_aero(id_aero), type, filename, iyr, klev_src, pt_ap, pt_b, pt_tmp, & |
| 213 |
|
✗ |
psurf_year(:,:,id_aero), load_year(:,:,id_aero)) |
| 214 |
|
✗ |
IF (.NOT. ALLOCATED(var_year)) THEN |
| 215 |
|
✗ |
ALLOCATE(var_year(klon, klev_src, 12, naero_spc), stat=ierr) |
| 216 |
|
✗ |
IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 5',1) |
| 217 |
|
|
ENDIF |
| 218 |
|
✗ |
var_year(:,:,:,id_aero) = pt_tmp(:,:,:) |
| 219 |
|
|
|
| 220 |
|
|
! Reading values corresponding to the preindustrial concentrations. |
| 221 |
|
✗ |
type='preind' |
| 222 |
|
|
CALL readaerosol(name_aero(id_aero), type, filename, iyr, pi_klev_src, pt_ap, pt_b, pt_tmp, & |
| 223 |
|
✗ |
pi_psurf_year(:,:,id_aero), pi_load_year(:,:,id_aero)) |
| 224 |
|
|
|
| 225 |
|
|
! klev_src must be the same in both files. |
| 226 |
|
|
! Also supposing pt_ap and pt_b to be the same in the 2 files without testing. |
| 227 |
|
✗ |
IF (pi_klev_src /= klev_src) THEN |
| 228 |
|
✗ |
WRITE(lunout,*) 'Error! All forcing files for the same aerosol must have the same vertical dimension' |
| 229 |
|
✗ |
WRITE(lunout,*) 'Aerosol : ', name_aero(id_aero) |
| 230 |
|
✗ |
CALL abort_physic('readaerosol_interp','Differnt vertical axes in aerosol forcing files',1) |
| 231 |
|
|
ENDIF |
| 232 |
|
|
|
| 233 |
|
✗ |
IF (.NOT. ALLOCATED(pi_var_year)) THEN |
| 234 |
|
✗ |
ALLOCATE(pi_var_year(klon, klev_src, 12, naero_spc), stat=ierr) |
| 235 |
|
✗ |
IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 6',1) |
| 236 |
|
|
ENDIF |
| 237 |
|
✗ |
pi_var_year(:,:,:,id_aero) = pt_tmp(:,:,:) |
| 238 |
|
|
|
| 239 |
|
✗ |
IF (debug) THEN |
| 240 |
|
✗ |
CALL writefield_phy('var_year_jan',var_year(:,:,1,id_aero),klev_src) |
| 241 |
|
✗ |
CALL writefield_phy('var_year_dec',var_year(:,:,12,id_aero),klev_src) |
| 242 |
|
✗ |
CALL writefield_phy('psurf_src',psurf_year(:,:,id_aero),1) |
| 243 |
|
✗ |
CALL writefield_phy('pi_psurf_src',pi_psurf_year(:,:,id_aero),1) |
| 244 |
|
✗ |
CALL writefield_phy('load_year_src',load_year(:,:,id_aero),1) |
| 245 |
|
✗ |
CALL writefield_phy('pi_load_year_src',pi_load_year(:,:,id_aero),1) |
| 246 |
|
|
ENDIF |
| 247 |
|
|
|
| 248 |
|
|
! Pointer no more useful, deallocate. |
| 249 |
|
✗ |
DEALLOCATE(pt_tmp) |
| 250 |
|
|
|
| 251 |
|
|
! Test if vertical interpolation will be needed. |
| 252 |
|
✗ |
IF (psurf_year(1,1,id_aero)==not_valid .OR. pi_psurf_year(1,1,id_aero)==not_valid ) THEN |
| 253 |
|
|
! Pressure=not_valid indicates old file format, see module readaerosol |
| 254 |
|
✗ |
vert_interp = .FALSE. |
| 255 |
|
|
|
| 256 |
|
|
! If old file format, both psurf_year and pi_psurf_year must be not_valid |
| 257 |
|
✗ |
IF ( psurf_year(1,1,id_aero) /= pi_psurf_year(1,1,id_aero) ) THEN |
| 258 |
|
✗ |
WRITE(lunout,*) 'Warning! All forcing files for the same aerosol must have the same structure' |
| 259 |
|
✗ |
CALL abort_physic('readaerosol_interp', 'The aerosol files have not the same format',1) |
| 260 |
|
|
ENDIF |
| 261 |
|
|
|
| 262 |
|
✗ |
IF (klev /= klev_src) THEN |
| 263 |
|
✗ |
WRITE(lunout,*) 'Old format of aerosol file do not allowed vertical interpolation' |
| 264 |
|
✗ |
CALL abort_physic('readaerosol_interp', 'Old aerosol file not possible',1) |
| 265 |
|
|
ENDIF |
| 266 |
|
|
|
| 267 |
|
|
ELSE |
| 268 |
|
✗ |
vert_interp = .TRUE. |
| 269 |
|
|
ENDIF |
| 270 |
|
|
|
| 271 |
|
|
! Calendar initialisation |
| 272 |
|
|
! |
| 273 |
|
✗ |
DO i = 2, 13 |
| 274 |
|
✗ |
month_len(i) = REAL(ioget_mon_len(year_cur, i-1)) |
| 275 |
|
✗ |
CALL ymds2ju(year_cur, i-1, 1, 0.0, month_start(i)) |
| 276 |
|
|
ENDDO |
| 277 |
|
✗ |
month_len(1) = REAL(ioget_mon_len(year_cur-1, 12)) |
| 278 |
|
✗ |
CALL ymds2ju(year_cur-1, 12, 1, 0.0, month_start(1)) |
| 279 |
|
✗ |
month_len(14) = REAL(ioget_mon_len(year_cur+1, 1)) |
| 280 |
|
✗ |
CALL ymds2ju(year_cur+1, 1, 1, 0.0, month_start(14)) |
| 281 |
|
✗ |
month_mid(:) = month_start (:) + month_len(:)/2. |
| 282 |
|
|
|
| 283 |
|
✗ |
if (debug) then |
| 284 |
|
✗ |
write(lunout,*)' month_len = ',month_len |
| 285 |
|
✗ |
write(lunout,*)' month_mid = ',month_mid |
| 286 |
|
|
endif |
| 287 |
|
|
|
| 288 |
|
|
ENDIF ! IF ( (first .OR. iday==0) .AND. lnewday ) THEN |
| 289 |
|
|
|
| 290 |
|
|
!**************************************************************************************** |
| 291 |
|
|
! - 2) Interpolate to the actual day. |
| 292 |
|
|
! - 3) Interpolate to the model vertical grid. |
| 293 |
|
|
! |
| 294 |
|
|
!**************************************************************************************** |
| 295 |
|
|
|
| 296 |
|
✗ |
IF (lnewday) THEN ! only if new day |
| 297 |
|
|
!**************************************************************************************** |
| 298 |
|
|
! 2) Interpolate to the actual day |
| 299 |
|
|
! |
| 300 |
|
|
!**************************************************************************************** |
| 301 |
|
|
! Find which months and days to use for time interpolation |
| 302 |
|
✗ |
nbr_tsteps = 12 |
| 303 |
|
|
IF (nbr_tsteps == 12) then |
| 304 |
|
✗ |
IF (jDay < month_mid(im+1)) THEN |
| 305 |
|
✗ |
im2=im-1 |
| 306 |
|
✗ |
day2 = month_mid(im2+1) |
| 307 |
|
✗ |
day1 = month_mid(im+1) |
| 308 |
|
✗ |
IF (im2 <= 0) THEN |
| 309 |
|
|
! the month is january, thus the month before december |
| 310 |
|
✗ |
im2=12 |
| 311 |
|
|
ENDIF |
| 312 |
|
|
ELSE |
| 313 |
|
|
! the second half of the month |
| 314 |
|
✗ |
im2=im+1 |
| 315 |
|
✗ |
day1 = month_mid(im+1) |
| 316 |
|
✗ |
day2 = month_mid(im2+1) |
| 317 |
|
✗ |
IF (im2 > 12) THEN |
| 318 |
|
|
! the month is december, the following thus january |
| 319 |
|
✗ |
im2=1 |
| 320 |
|
|
ENDIF |
| 321 |
|
|
ENDIF |
| 322 |
|
|
ELSE IF (nbr_tsteps == 14) then |
| 323 |
|
|
im = im + 1 |
| 324 |
|
|
IF (jDay < month_mid(im)) THEN |
| 325 |
|
|
! in the first half of the month use month before and actual month |
| 326 |
|
|
im2=im-1 |
| 327 |
|
|
day2 = month_mid(im2) |
| 328 |
|
|
day1 = month_mid(im) |
| 329 |
|
|
ELSE |
| 330 |
|
|
! the second half of the month |
| 331 |
|
|
im2=im+1 |
| 332 |
|
|
day1 = month_mid(im) |
| 333 |
|
|
day2 = month_mid(im2) |
| 334 |
|
|
ENDIF |
| 335 |
|
|
ELSE |
| 336 |
|
|
CALL abort_physic('readaerosol_interp', 'number of months undefined',1) |
| 337 |
|
|
ENDIF |
| 338 |
|
✗ |
if (debug) then |
| 339 |
|
✗ |
write(lunout,*)' jDay, day1, day2, im, im2 = ', jDay, day1, day2, im, im2 |
| 340 |
|
|
endif |
| 341 |
|
|
|
| 342 |
|
|
|
| 343 |
|
|
! Time interpolation, still on vertical source grid |
| 344 |
|
✗ |
ALLOCATE(tmp1(klon,klev_src), tmp2(klon,klev_src),stat=ierr) |
| 345 |
|
✗ |
IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 7',1) |
| 346 |
|
|
|
| 347 |
|
✗ |
ALLOCATE(pplay_src(klon,klev_src), stat=ierr) |
| 348 |
|
✗ |
IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 8',1) |
| 349 |
|
|
|
| 350 |
|
|
|
| 351 |
|
✗ |
DO k=1,klev_src |
| 352 |
|
✗ |
DO i=1,klon |
| 353 |
|
|
tmp1(i,k) = & |
| 354 |
|
|
var_year(i,k,im2,id_aero) - (jDay-day2)/(day1-day2) * & |
| 355 |
|
✗ |
(var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero)) |
| 356 |
|
|
|
| 357 |
|
|
tmp2(i,k) = & |
| 358 |
|
|
pi_var_year(i,k,im2,id_aero) - (jDay-day2)/(day1-day2) * & |
| 359 |
|
✗ |
(pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero)) |
| 360 |
|
|
ENDDO |
| 361 |
|
|
ENDDO |
| 362 |
|
|
|
| 363 |
|
|
! Time interpolation for pressure at surface, still on vertical source grid |
| 364 |
|
✗ |
DO i=1,klon |
| 365 |
|
|
psurf_day(i) = & |
| 366 |
|
|
psurf_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * & |
| 367 |
|
✗ |
(psurf_year(i,im2,id_aero) - psurf_year(i,im,id_aero)) |
| 368 |
|
|
|
| 369 |
|
|
pi_psurf_day(i) = & |
| 370 |
|
|
pi_psurf_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * & |
| 371 |
|
✗ |
(pi_psurf_year(i,im2,id_aero) - pi_psurf_year(i,im,id_aero)) |
| 372 |
|
|
ENDDO |
| 373 |
|
|
|
| 374 |
|
|
! Time interpolation for the load, still on vertical source grid |
| 375 |
|
✗ |
DO i=1,klon |
| 376 |
|
|
load_src(i) = & |
| 377 |
|
|
load_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * & |
| 378 |
|
✗ |
(load_year(i,im2,id_aero) - load_year(i,im,id_aero)) |
| 379 |
|
|
|
| 380 |
|
|
pi_load_src(i) = & |
| 381 |
|
|
pi_load_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * & |
| 382 |
|
✗ |
(pi_load_year(i,im2,id_aero) - pi_load_year(i,im,id_aero)) |
| 383 |
|
|
ENDDO |
| 384 |
|
|
|
| 385 |
|
|
!**************************************************************************************** |
| 386 |
|
|
! 3) Interpolate to the model vertical grid (target grid) |
| 387 |
|
|
! |
| 388 |
|
|
!**************************************************************************************** |
| 389 |
|
|
|
| 390 |
|
✗ |
IF (vert_interp) THEN |
| 391 |
|
|
|
| 392 |
|
|
! - Interpolate variable tmp1 (on source grid) to var_day (on target grid) |
| 393 |
|
|
!******************************************************************************** |
| 394 |
|
|
! a) calculate pression at vertical levels for the source grid using the |
| 395 |
|
|
! hybrid-sigma coordinates ap and b and the surface pressure, variables from file. |
| 396 |
|
✗ |
DO k = 1, klev_src |
| 397 |
|
✗ |
DO i = 1, klon |
| 398 |
|
✗ |
pplay_src(i,k)= pt_ap(k) + pt_b(k)*psurf_day(i) |
| 399 |
|
|
ENDDO |
| 400 |
|
|
ENDDO |
| 401 |
|
|
|
| 402 |
|
✗ |
IF (debug) THEN |
| 403 |
|
✗ |
CALL writefield_phy('psurf_day_src',psurf_day(:),1) |
| 404 |
|
✗ |
CALL writefield_phy('pplay_src',pplay_src(:,:),klev_src) |
| 405 |
|
✗ |
CALL writefield_phy('pplay',pplay(:,:),klev) |
| 406 |
|
✗ |
CALL writefield_phy('day_src',tmp1,klev_src) |
| 407 |
|
✗ |
CALL writefield_phy('pi_day_src',tmp2,klev_src) |
| 408 |
|
|
ENDIF |
| 409 |
|
|
|
| 410 |
|
|
! b) vertical interpolation on pressure leveles |
| 411 |
|
|
CALL pres2lev(tmp1(:,:), var_day(:,:,id_aero), klev_src, klev, pplay_src, pplay, & |
| 412 |
|
✗ |
1, klon, .FALSE.) |
| 413 |
|
|
|
| 414 |
|
✗ |
IF (debug) CALL writefield_phy('day_tgt',var_day(:,:,id_aero),klev) |
| 415 |
|
|
|
| 416 |
|
|
! c) adjust to conserve total aerosol mass load in the vertical pillar |
| 417 |
|
|
! Calculate the load in the actual pillar and compare with the load |
| 418 |
|
|
! read from aerosol file. |
| 419 |
|
|
|
| 420 |
|
|
! Find the pressure difference in each model layer |
| 421 |
|
✗ |
DO k = 1, klev |
| 422 |
|
✗ |
DO i = 1, klon |
| 423 |
|
✗ |
delp(i,k) = paprs(i,k) - paprs (i,k+1) |
| 424 |
|
|
ENDDO |
| 425 |
|
|
ENDDO |
| 426 |
|
|
|
| 427 |
|
|
! Find the mass load in the actual pillar, on target grid |
| 428 |
|
✗ |
load_tgt(:) = 0. |
| 429 |
|
✗ |
DO k= 1, klev |
| 430 |
|
✗ |
DO i = 1, klon |
| 431 |
|
✗ |
zrho = pplay(i,k)/t_seri(i,k)/RD ! [kg/m3] |
| 432 |
|
✗ |
volm = var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg] |
| 433 |
|
✗ |
load_tgt(i) = load_tgt(i) + volm *delp(i,k)/RG |
| 434 |
|
|
ENDDO |
| 435 |
|
|
ENDDO |
| 436 |
|
|
|
| 437 |
|
|
! Adjust, uniform |
| 438 |
|
✗ |
DO k = 1, klev |
| 439 |
|
✗ |
DO i = 1, klon |
| 440 |
|
✗ |
var_day(i,k,id_aero) = var_day(i,k,id_aero)*load_src(i)/max(1.e-30,load_tgt(i)) |
| 441 |
|
|
ENDDO |
| 442 |
|
|
ENDDO |
| 443 |
|
|
|
| 444 |
|
✗ |
IF (debug) THEN |
| 445 |
|
✗ |
load_tgt_test(:) = 0. |
| 446 |
|
✗ |
DO k= 1, klev |
| 447 |
|
✗ |
DO i = 1, klon |
| 448 |
|
✗ |
zrho = pplay(i,k)/t_seri(i,k)/RD ! [kg/m3] |
| 449 |
|
✗ |
volm = var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg] |
| 450 |
|
✗ |
load_tgt_test(i) = load_tgt_test(i) + volm*delp(i,k)/RG |
| 451 |
|
|
ENDDO |
| 452 |
|
|
ENDDO |
| 453 |
|
|
|
| 454 |
|
✗ |
CALL writefield_phy('day_tgt2',var_day(:,:,id_aero),klev) |
| 455 |
|
✗ |
CALL writefield_phy('load_tgt',load_tgt(:),1) |
| 456 |
|
✗ |
CALL writefield_phy('load_tgt_test',load_tgt_test(:),1) |
| 457 |
|
✗ |
CALL writefield_phy('load_src',load_src(:),1) |
| 458 |
|
|
ENDIF |
| 459 |
|
|
|
| 460 |
|
|
! - Interpolate variable tmp2 (source grid) to pi_var_day (target grid) |
| 461 |
|
|
!******************************************************************************** |
| 462 |
|
|
! a) calculate pression at vertical levels at source grid |
| 463 |
|
✗ |
DO k = 1, klev_src |
| 464 |
|
✗ |
DO i = 1, klon |
| 465 |
|
✗ |
pplay_src(i,k)= pt_ap(k) + pt_b(k)*pi_psurf_day(i) |
| 466 |
|
|
ENDDO |
| 467 |
|
|
ENDDO |
| 468 |
|
|
|
| 469 |
|
✗ |
IF (debug) THEN |
| 470 |
|
✗ |
CALL writefield_phy('pi_psurf_day_src',pi_psurf_day(:),1) |
| 471 |
|
✗ |
CALL writefield_phy('pi_pplay_src',pplay_src(:,:),klev_src) |
| 472 |
|
|
ENDIF |
| 473 |
|
|
|
| 474 |
|
|
! b) vertical interpolation on pressure leveles |
| 475 |
|
|
CALL pres2lev(tmp2(:,:), pi_var_day(:,:,id_aero), klev_src, klev, pplay_src, pplay, & |
| 476 |
|
✗ |
1, klon, .FALSE.) |
| 477 |
|
|
|
| 478 |
|
✗ |
IF (debug) CALL writefield_phy('pi_day_tgt',pi_var_day(:,:,id_aero),klev) |
| 479 |
|
|
|
| 480 |
|
|
! c) adjust to conserve total aerosol mass load in the vertical pillar |
| 481 |
|
|
! Calculate the load in the actual pillar and compare with the load |
| 482 |
|
|
! read from aerosol file. |
| 483 |
|
|
|
| 484 |
|
|
! Find the load in the actual pillar, on target grid |
| 485 |
|
✗ |
load_tgt(:) = 0. |
| 486 |
|
✗ |
DO k = 1, klev |
| 487 |
|
✗ |
DO i = 1, klon |
| 488 |
|
✗ |
zrho = pplay(i,k)/t_seri(i,k)/RD ! [kg/m3] |
| 489 |
|
✗ |
volm = pi_var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg] |
| 490 |
|
✗ |
load_tgt(i) = load_tgt(i) + volm*delp(i,k)/RG |
| 491 |
|
|
ENDDO |
| 492 |
|
|
ENDDO |
| 493 |
|
|
|
| 494 |
|
✗ |
DO k = 1, klev |
| 495 |
|
✗ |
DO i = 1, klon |
| 496 |
|
✗ |
pi_var_day(i,k,id_aero) = pi_var_day(i,k,id_aero)*pi_load_src(i)/max(1.e-30,load_tgt(i)) |
| 497 |
|
|
ENDDO |
| 498 |
|
|
ENDDO |
| 499 |
|
|
|
| 500 |
|
✗ |
IF (debug) THEN |
| 501 |
|
✗ |
load_tgt_test(:) = 0. |
| 502 |
|
✗ |
DO k = 1, klev |
| 503 |
|
✗ |
DO i = 1, klon |
| 504 |
|
✗ |
zrho = pplay(i,k)/t_seri(i,k)/RD ! [kg/m3] |
| 505 |
|
✗ |
volm = pi_var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg] |
| 506 |
|
✗ |
load_tgt_test(i) = load_tgt_test(i) + volm*delp(i,k)/RG |
| 507 |
|
|
ENDDO |
| 508 |
|
|
ENDDO |
| 509 |
|
✗ |
CALL writefield_phy('pi_day_tgt2',pi_var_day(:,:,id_aero),klev) |
| 510 |
|
✗ |
CALL writefield_phy('pi_load_tgt',load_tgt(:),1) |
| 511 |
|
✗ |
CALL writefield_phy('pi_load_tgt_test',load_tgt_test(:),1) |
| 512 |
|
✗ |
CALL writefield_phy('pi_load_src',pi_load_src(:),1) |
| 513 |
|
|
ENDIF |
| 514 |
|
|
|
| 515 |
|
|
|
| 516 |
|
|
ELSE ! No vertical interpolation done |
| 517 |
|
|
|
| 518 |
|
✗ |
var_day(:,:,id_aero) = tmp1(:,:) |
| 519 |
|
✗ |
pi_var_day(:,:,id_aero) = tmp2(:,:) |
| 520 |
|
|
|
| 521 |
|
|
ENDIF ! vert_interp |
| 522 |
|
|
|
| 523 |
|
|
|
| 524 |
|
|
! Deallocation |
| 525 |
|
✗ |
DEALLOCATE(tmp1, tmp2, pplay_src, stat=ierr) |
| 526 |
|
|
|
| 527 |
|
|
!**************************************************************************************** |
| 528 |
|
|
! 4) Test for negative mass values |
| 529 |
|
|
! |
| 530 |
|
|
!**************************************************************************************** |
| 531 |
|
✗ |
IF (MINVAL(var_day(:,:,id_aero)) < 0.) THEN |
| 532 |
|
✗ |
DO k=1,klev |
| 533 |
|
✗ |
DO i=1,klon |
| 534 |
|
|
! Test for var_day |
| 535 |
|
✗ |
IF (var_day(i,k,id_aero) < 0.) THEN |
| 536 |
|
✗ |
IF (jDay-day2 < 0.) WRITE(lunout,*) 'jDay-day2=',jDay-day2 |
| 537 |
|
✗ |
IF (var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero) < 0.) THEN |
| 538 |
|
✗ |
WRITE(lunout,*) trim(name_aero(id_aero)),'(i,k,im2)-', & |
| 539 |
|
✗ |
trim(name_aero(id_aero)),'(i,k,im)=', & |
| 540 |
|
✗ |
var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero) |
| 541 |
|
|
ENDIF |
| 542 |
|
✗ |
WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero) |
| 543 |
|
✗ |
WRITE(lunout,*) 'day1, day2, jDay = ', day1, day2, jDay |
| 544 |
|
✗ |
CALL abort_physic('readaerosol_interp','Error in interpolation 1',1) |
| 545 |
|
|
ENDIF |
| 546 |
|
|
ENDDO |
| 547 |
|
|
ENDDO |
| 548 |
|
|
ENDIF |
| 549 |
|
|
|
| 550 |
|
✗ |
IF (MINVAL(pi_var_day(:,:,id_aero)) < 0. ) THEN |
| 551 |
|
✗ |
DO k=1, klev |
| 552 |
|
✗ |
DO i=1,klon |
| 553 |
|
|
! Test for pi_var_day |
| 554 |
|
✗ |
IF (pi_var_day(i,k,id_aero) < 0.) THEN |
| 555 |
|
✗ |
IF (jDay-day2 < 0.) WRITE(lunout,*) 'jDay-day2=',jDay-day2 |
| 556 |
|
✗ |
IF (pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero) < 0.) THEN |
| 557 |
|
✗ |
WRITE(lunout,*) trim(name_aero(id_aero)),'(i,k,im2)-', & |
| 558 |
|
✗ |
trim(name_aero(id_aero)),'(i,k,im)=', & |
| 559 |
|
✗ |
pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero) |
| 560 |
|
|
ENDIF |
| 561 |
|
|
|
| 562 |
|
✗ |
WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero) |
| 563 |
|
✗ |
CALL abort_physic('readaerosol_interp','Error in interpolation 2',1) |
| 564 |
|
|
ENDIF |
| 565 |
|
|
ENDDO |
| 566 |
|
|
ENDDO |
| 567 |
|
|
ENDIF |
| 568 |
|
|
|
| 569 |
|
|
ENDIF ! lnewday |
| 570 |
|
|
|
| 571 |
|
|
!**************************************************************************************** |
| 572 |
|
|
! Copy output from saved variables |
| 573 |
|
|
! |
| 574 |
|
|
!**************************************************************************************** |
| 575 |
|
|
|
| 576 |
|
✗ |
mass_out(:,:) = var_day(:,:,id_aero) |
| 577 |
|
✗ |
pi_mass_out(:,:) = pi_var_day(:,:,id_aero) |
| 578 |
|
|
|
| 579 |
|
✗ |
END SUBROUTINE readaerosol_interp |
| 580 |
|
|
|