! ! $Header$ ! SUBROUTINE readaerosol (id_aero, r_day, first, massvar_p) USE dimphy, ONLY : klev USE mod_grid_phy_lmdz, klon=>klon_glo USE mod_phys_lmdz_para IMPLICIT NONE ! Content: ! -------- ! This routine reads in monthly mean values of massvar aerosols and ! returns a linearly interpolated dayly-mean field. ! ! ! Author: ! ------- ! Johannes Quaas (quaas@lmd.jussieu.fr) ! Celine Deandreis & Anne Cozic ! 19/02/09 ! ! ! Include-files: ! -------------- INCLUDE "YOMCST.h" INCLUDE "chem.h" INCLUDE "dimensions.h" INCLUDE "temps.h" INCLUDE "clesphys.h" INCLUDE "iniprint.h" ! ! Input: ! ------ INTEGER, INTENT(in) :: id_aero REAL, INTENT(in) :: r_day ! Day of integration LOGICAL, INTENT(in) :: first ! First timestep ! (and therefore initialization necessary) ! ! Output: ! ------- REAL, INTENT(out) :: massvar_p(klon_omp,klev) ! Mass of massvar (monthly mean data,from file) [ug AIBCM/m3] ! ! Local Variables: ! ---------------- INTEGER :: i, ig, k, it INTEGER :: j, iday, iyr, iyr1, iyr2 INTEGER :: im, day1, day2, im2 INTEGER, PARAMETER :: ny=jjm+1 CHARACTER(len=4) :: cyear CHARACTER(len=7),DIMENSION(8) :: name_aero REAL :: var_1(iim, jjm+1, klev, 12) REAL, DIMENSION(klon,klev) :: massvar REAL, DIMENSION(iim, jjm+1, klev, 12) :: var_2 ! The massvar distributions REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: var ! VAR in right dimension REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: var_out !$OMP THREADPRIVATE(var,var_out) LOGICAL :: lnewday LOGICAL, PARAMETER :: lonlyone=.FALSE. LOGICAL,SAVE :: first2=.TRUE. !$OMP THREADPRIVATE(first2) ! Variable pour definir a partir d'un numero d'aerosol, son nom name_aero(1) = "SSSSM " name_aero(2) = "ASSSM " name_aero(3) = "ASBCM " name_aero(4) = "ASPOMM " name_aero(5) = "SO4 " name_aero(6) = "CIDUSTM" name_aero(7) = "AIBCM " name_aero(8) = "AIPOMM " !$OMP MASTER IF (first2) THEN ALLOCATE( var(klon, klev, 12,8) ) ALLOCATE( var_out(klon, klev,8)) first2=.FALSE. IF (aer_type /= 'actuel ' .AND. aer_type /= 'preind ' .AND. & aer_type /= 'scenario') THEN WRITE(lunout,*)' *** Warning ***' WRITE(lunout,*)'Option aer_type non prevu pour les aerosols = ', & aer_type CALL abort_gcm('readaerosol','Error : aer_type parameter not accepted',1) ENDIF ENDIF IF (is_mpi_root) THEN iday = INT(r_day) ! Get the year of the run iyr = iday/360 ! Get the day of the actual year: iday = iday-iyr*360 ! 0.02 is about 0.5/24, namly less than half an hour lnewday = (r_day-FLOAT(iday).LT.0.02) ! --------------------------------------------- ! All has to be done only, if a new day begins! ! --------------------------------------------- IF (lnewday.OR.first) THEN im = iday/30 +1 ! the actual month ! annee_ref is the initial year (defined in temps.h) iyr = iyr + annee_ref ! Do I have to read new data? (Is this the first day of a year?) IF (first.OR.iday.EQ.1.) THEN ! Initialize values DO it=1,12 DO k=1,klev DO i=1,klon var(i,k,it,id_aero)=0. ENDDO ENDDO ENDDO IF (aer_type == 'actuel ') then cyear='1980' CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero)) ELSE IF (aer_type == 'preind ') THEN cyear='.nat' CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero)) ELSE ! aer_type=scenario IF (iyr .LT. 1850) THEN cyear='.nat' WRITE(lunout,*) 'get_aero iyr=', iyr,' ',cyear CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero)) ELSE IF (iyr .GE. 2100) THEN cyear='2100' WRITE(lunout,*) 'get_aero iyr=', iyr,' ',cyear CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero)) ELSE ! Read data from 2 decades ! a) from actual 10-yr-period IF (iyr.LT.1900) THEN iyr1 = 1850 iyr2 = 1900 ELSE IF (iyr.GE.1900.AND.iyr.LT.1920) THEN iyr1 = 1900 iyr2 = 1920 ELSE iyr1 = INT(iyr/10)*10 iyr2 = INT(1+iyr/10)*10 ENDIF WRITE(cyear,'(I4)') iyr1 WRITE(lunout,*) 'get_aero iyr=', iyr,' ',cyear CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero)) ! If to read two decades: IF (.NOT.lonlyone) THEN ! b) from the next following one WRITE(cyear,'(I4)') iyr2 WRITE(lunout,*) 'get_aero iyr=', iyr,' ',cyear CALL get_aero_fromfile(cyear, var_2, name_aero(id_aero)) ! Interpolate linarily to the actual year: DO it=1,12 DO k=1,klev DO j=1,jjm DO i=1,iim var_1(i,j,k,it) = & var_1(i,j,k,it) - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1) * & (var_1(i,j,k,it) - var_2(i,j,k,it)) ENDDO ENDDO ENDDO ENDDO ENDIF ! lonlyone ENDIF ! iyr .LT. 1850 ENDIF ! aer_type ! Transform the horizontal 2D-field into the physics-field ! (Also the levels and the latitudes have to be inversed) DO it=1,12 DO k=1, klev ! a) at the poles, use the zonal mean: DO i=1,iim ! North pole var(1,k,it,id_aero) = & var(1,k,it,id_aero)+ var_1(i,jjm+1,klev+1-k,it) ! South pole var(klon,k,it,id_aero)= & var(klon,k,it,id_aero)+ var_1(i,1,klev+1-k,it) ENDDO var(1,k,it,id_aero) = var(1,k,it,id_aero)/FLOAT(iim) var(klon,k,it,id_aero)= var(klon,k,it,id_aero)/FLOAT(iim) ! b) the values between the poles: ig=1 DO j=2,jjm DO i=1,iim ig=ig+1 IF (ig.GT.klon) THEN WRITE(lunout,*) 'Attention dans readaerosol', & name_aero(id_aero), 'ig > klon', ig, klon ENDIF var(ig,k,it,id_aero) = var_1(i,jjm+1+1-j,klev+1-k,it) ENDDO ENDDO IF (ig.NE.klon-1) THEN PRINT *,'Error in readaerosol', name_aero(id_aero), 'conversion' CALL abort_gcm('readaerosol','Error in readaerosol 1',1) ENDIF ENDDO ! Loop over k (vertical) ENDDO ! Loop over it (months) ENDIF ! Had to read new data? ! Interpolate to actual day: IF (iday.LT.im*30-15) THEN ! in the first half of the month use month before and actual month im2=im-1 day2 = im2*30-15 day1 = im2*30+15 IF (im2.LE.0) THEN ! the month is january, thus the month before december im2=12 ENDIF DO k=1,klev DO i=1,klon massvar(i,k) = & var(i,k,im2,id_aero)- FLOAT(iday-day2)/FLOAT(day1-day2) * & (var(i,k,im2,id_aero) - var(i,k,im,id_aero)) IF (massvar(i,k).LT.0.) THEN IF (iday-day2.LT.0.) WRITE(lunout,*) 'iday-day2',iday-day2 IF (var(i,k,im2,id_aero) - var(i,k,im,id_aero).LT.0.) THEN PRINT*, name_aero(id_aero),'(i,k,im2)- ', & name_aero(id_aero),'(i,k,im)', & var(i,k,im2,id_aero) - var(i,k,im,id_aero) ENDIF IF (day1-day2.LT.0.) WRITE(lunout,*) 'day1-day2',day1-day2 WRITE(lunout,*) 'stop',name_aero(id_aero) CALL abort_gcm('readaerosol','Error in interpolation 2',1) ENDIF ENDDO ENDDO ELSE ! the second half of the month im2=im+1 IF (im2.GT.12) THEN ! the month is december, the following thus january im2=1 ENDIF day2 = im*30-15 day1 = im*30+15 DO k=1,klev DO i=1,klon massvar(i,k) = & var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * & (var(i,k,im2,id_aero) - var(i,k,im,id_aero)) IF (massvar(i,k).LT.0.) THEN IF (iday-day2.LT.0.) & WRITE(lunout,*) 'iday-day2',iday-day2 IF (var(i,k,im2,id_aero) - var(i,k,im,id_aero).LT.0.) THEN WRITE(lunout,*) name_aero(id_aero),'(i,k,im2)-', & name_aero(id_aero),'(i,k,im)', & var(i,k,im2,id_aero) - var(i,k,im,id_aero) ENDIF IF (day1-day2.LT.0.) & WRITE(lunout,*) 'day1-day2',day1-day2 WRITE(lunout,*) 'stop',name_aero(id_aero) CALL abort_gcm('readaerosol','Error in interpolation 3',1) ENDIF ENDDO ENDDO ENDIF ! The massvar concentration [molec cm-3] is read in. ! Convert it into mass [ug VAR/m3] ! masse_var in [g/mol], n_avogadro in [molec/mol] ! The sulfate mass [ug VAR/m3] is read in. DO k=1,klev DO i=1,klon var_out(i,k,id_aero) = massvar(i,k) IF (var_out(i,k,id_aero).LT.0) THEN PRINT*, 'Attention il y a un probleme dans readaerosol', & name_aero(id_aero), 'la masse est negative' CALL abort_gcm('readaerosol','Error in readaerosol 4',1) ENDIF ENDDO ENDDO ELSE ! if no new day, use old data: DO k=1,klev DO i=1,klon massvar(i,k) = var_out(i,k,id_aero) IF (var_out(i,k,id_aero).LT.0) THEN PRINT*, 'Attention il y a un probleme dans readaerosol', & name_aero(id_aero), 'la masse est negative' CALL abort_gcm('readaerosol','Error in readaerosol 5',1) ENDIF ENDDO ENDDO ENDIF ! Did I have to do anything (was it a new day?) ENDIF ! phy_rank==0 !$OMP END MASTER CALL Scatter(massvar,massvar_p) END SUBROUTINE readaerosol !----------------------------------------------------------------------------- ! Read in /calculate pre-industrial values of sulfate !----------------------------------------------------------------------------- SUBROUTINE readaerosol_preind (id_aero, r_day, first, pi_massvar_p) USE dimphy, ONLY : klev USE mod_grid_phy_lmdz, klon=>klon_glo USE mod_phys_lmdz_para IMPLICIT NONE ! Content: ! -------- ! This routine reads in monthly mean values of massvar aerosols and ! returns a linearly interpolated dayly-mean field. ! ! It does so for the preindustriel values of the massvar, to a large part ! analogous to the routine readmassvar above. ! ! Only Pb: Variables must be saved and don t have to be overwritten! ! ! Author: ! ------- ! Celine Deandreis & Anne Cozic (LSCE) 2009 ! Johannes Quaas (quaas@lmd.jussieu.fr) ! 26/06/01 ! ! Modifications: ! -------------- ! see above ! Include-files: ! -------------- INCLUDE "YOMCST.h" INCLUDE "chem.h" INCLUDE "dimensions.h" INCLUDE "temps.h" INCLUDE "iniprint.h" ! Input: ! ------ INTEGER, INTENT(in) :: id_aero REAL, INTENT(in) :: r_day ! Day of integration LOGICAL, INTENT(in) :: first ! First timestep (and therefore initialization necessary) ! ! Output: ! ------- REAL, DIMENSION(klon_omp,klev), INTENT(out) :: pi_massvar_p ! Number conc. massvar (monthly mean data, from file) ! ! Local Variables: ! ---------------- INTEGER :: i, ig, k, it INTEGER :: j, iday, iyr INTEGER, PARAMETER :: ny=jjm+1 INTEGER :: im, day1, day2, im2 REAL, DIMENSION(iim, jjm+1, klev, 12) :: pi_var_1 REAL, DIMENSION(klon,klev) :: pi_massvar REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: pi_var ! VAR in right dimension REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: pi_var_out !$OMP THREADPRIVATE(pi_var,pi_var_out) CHARACTER(len=4) :: cyear CHARACTER(len=7),DIMENSION(8) :: name_aero LOGICAL :: lnewday LOGICAL,SAVE :: first2=.TRUE. !$OMP THREADPRIVATE(first2) ! Variable pour definir a partir d'un numero d'aerosol, son nom name_aero(1) = "SSSSM " name_aero(2) = "ASSSM " name_aero(3) = "ASBCM " name_aero(4) = "ASPOMM " name_aero(5) = "SO4 " name_aero(6) = "CIDUSTM" name_aero(7) = "AIBCM " name_aero(8) = "AIPOMM " !$OMP MASTER IF (first2) THEN ALLOCATE( pi_var(klon, klev, 12,8) ) ALLOCATE( pi_var_out(klon, klev,8)) first2=.FALSE. ENDIF IF (is_mpi_root) THEN iday = INT(r_day) ! Get the year of the run iyr = iday/360 ! Get the day of the actual year: iday = iday-iyr*360 ! 0.02 is about 0.5/24, namly less than half an hour lnewday = (r_day-FLOAT(iday).LT.0.02) ! --------------------------------------------- ! All has to be done only, if a new day begins! ! --------------------------------------------- IF (lnewday.OR.first) THEN im = iday/30 +1 ! the actual month ! annee_ref is the initial year (defined in temps.h) iyr = iyr + annee_ref IF (first) THEN cyear='.nat' CALL get_aero_fromfile(cyear,pi_var_1, name_aero(id_aero)) ! Transform the horizontal 2D-field into the physics-field ! (Also the levels and the latitudes have to be inversed) ! Initialize field DO it=1,12 DO k=1,klev DO i=1,klon pi_var(i,k,it,id_aero)=0. ENDDO ENDDO ENDDO WRITE(lunout,*) 'preind: finished reading', FLOAT(iim) DO it=1,12 DO k=1, klev ! a) at the poles, use the zonal mean: DO i=1,iim ! North pole pi_var(1,k,it,id_aero) = & pi_var(1,k,it,id_aero) + pi_var_1(i,jjm+1,klev+1-k,it) ! South pole pi_var(klon,k,it,id_aero) = & pi_var(klon,k,it,id_aero) + pi_var_1(i,1,klev+1-k,it) ENDDO pi_var(1,k,it,id_aero) = pi_var(1,k,it,id_aero)/FLOAT(iim) pi_var(klon,k,it,id_aero) = pi_var(klon,k,it,id_aero)/FLOAT(iim) ! b) the values between the poles: ig=1 DO j=2,jjm DO i=1,iim ig=ig+1 IF (ig.GT.klon) THEN WRITE(lunout,*) 'Attention dans readaerosol_preind', & name_aero(id_aero), 'ig > klon', ig, klon ENDIF pi_var(ig,k,it,id_aero) = & pi_var_1(i,jjm+1+1-j,klev+1-k,it) ENDDO ENDDO IF (ig.NE.klon-1) THEN WRITE(lunout,*) 'Error in readaerosol_preind (', name_aero(id_aero), ' conversion)' CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 1',1) ENDIF ENDDO ! Loop over k (vertical) ENDDO ! Loop over it (months) ENDIF ! Had to read new data? ! Interpolate to actual day: IF (iday.LT.im*30-15) THEN ! in the first half of the month use month before and actual month im2=im-1 day1 = im2*30+15 day2 = im2*30-15 IF (im2.LE.0) THEN ! the month is january, thus the month before december im2=12 ENDIF DO k=1,klev DO i=1,klon pi_massvar(i,k) = & pi_var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * & (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero)) IF (pi_massvar(i,k).LT.0.) THEN IF (iday-day2.LT.0.) WRITE(lunout,*) 'iday-day2',iday-day2 IF (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero).LT.0.) THEN WRITE(lunout,*) 'pi_',name_aero(id_aero),'(i,k,im2) - pi_', & name_aero(id_aero),'(i,k,im)', & pi_var(i,k,im2,id_aero) -pi_var(i,k,im,id_aero) ENDIF IF (day1-day2.LT.0.)WRITE(lunout,*) 'day1-day2',day1-day2 PRINT *, 'pi_',name_aero(id_aero) CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 2',1) ENDIF ENDDO ENDDO ELSE ! the second half of the month im2=im+1 day1 = im*30+15 IF (im2.GT.12) THEN ! the month is december, the following thus january im2=1 ENDIF day2 = im*30-15 DO k=1,klev DO i=1,klon pi_massvar(i,k) = & pi_var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * & (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero)) IF (pi_massvar(i,k).LT.0.) THEN IF (iday-day2.LT.0.) WRITE(lunout,*) 'iday-day2',iday-day2 IF (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero).LT.0.) THEN WRITE(lunout,*) 'pi_', name_aero(id_aero),'(i,k,im2) - pi_',& name_aero(id_aero), '(i,k,im)',& pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero) ENDIF IF (day1-day2.LT.0.) & WRITE(lunout,*) 'day1-day2',day1-day2 PRINT *, 'pi_',name_aero(id_aero) CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 3',1) ENDIF ENDDO ENDDO ENDIF ! The massvar concentration [molec cm-3] is read in. ! Convert it into mass [ug VAR/m3] ! masse_var in [g/mol], n_avogadro in [molec/mol] DO k=1,klev DO i=1,klon pi_var_out(i,k,id_aero) = pi_massvar(i,k) ENDDO ENDDO ELSE ! If no new day, use old data: DO k=1,klev DO i=1,klon pi_massvar(i,k) = pi_var_out(i,k,id_aero) ENDDO ENDDO ENDIF ! Was this the beginning of a new day? ENDIF ! is_mpi_root !$OMP END MASTER CALL Scatter(pi_massvar,pi_massvar_p) END SUBROUTINE readaerosol_preind !----------------------------------------------------------------------------- ! Routine for reading VAR data from files !----------------------------------------------------------------------------- SUBROUTINE get_aero_fromfile (cyr, var, varname) USE dimphy IMPLICIT NONE INCLUDE "netcdf.inc" INCLUDE "dimensions.h" INCLUDE "iniprint.h" CHARACTER(len=4), INTENT(in) :: cyr REAL, DIMENSION(iim, jjm+1, klev, 12), INTENT(out) :: var CHARACTER(len=*), INTENT(in) :: varname CHARACTER(len=30) :: fname CHARACTER(len=30) :: cvar INTEGER, DIMENSION(3) :: START, COUNT INTEGER :: STATUS, NCID, VARID INTEGER :: imth, i, j, k INTEGER, PARAMETER :: ny=jjm+1 REAL, DIMENSION(iim, ny, klev) :: varmth ! ! ! fname = TRIM(varname)//'.run'//cyr//'.cdf' WRITE(lunout,*) 'reading ', TRIM(fname) STATUS = NF_OPEN (TRIM(fname), NF_NOWRITE, NCID) IF (STATUS .NE. NF_NOERR) THEN WRITE(lunout,*) 'err in open ',status CALL abort_gcm('get_aero_fromfile','Error in opening file',1) ENDIF DO imth=1, 12 IF (imth.EQ.1) THEN cvar=TRIM(varname)//'JAN' ELSEIF (imth.EQ.2) THEN cvar=TRIM(varname)//'FEB' ELSEIF (imth.EQ.3) THEN cvar=TRIM(varname)//'MAR' ELSEIF (imth.EQ.4) THEN cvar=TRIM(varname)//'APR' ELSEIF (imth.EQ.5) THEN cvar=TRIM(varname)//'MAY' ELSEIF (imth.EQ.6) THEN cvar=TRIM(varname)//'JUN' ELSEIF (imth.EQ.7) THEN cvar=TRIM(varname)//'JUL' ELSEIF (imth.EQ.8) THEN cvar=TRIM(varname)//'AUG' ELSEIF (imth.EQ.9) THEN cvar=TRIM(varname)//'SEP' ELSEIF (imth.EQ.10) THEN cvar=TRIM(varname)//'OCT' ELSEIF (imth.EQ.11) THEN cvar=TRIM(varname)//'NOV' ELSEIF (imth.EQ.12) THEN cvar=TRIM(varname)//'DEC' ENDIF start(1)=1 start(2)=1 start(3)=1 COUNT(1)=iim COUNT(2)=ny COUNT(3)=klev STATUS = NF_INQ_VARID (NCID, TRIM(cvar), VARID) WRITE(lunout,*) ncid,imth,TRIM(cvar), varid IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in read ',status #ifdef NC_DOUBLE status = NF_GET_VARA_DOUBLE(NCID, VARID, START, COUNT,varmth) #else status = NF_GET_VARA_REAL(NCID, VARID, START, COUNT, varmth) #endif IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in read data',status DO k=1,klev DO j=1,jjm+1 DO i=1,iim IF (varmth(i,j,k).LT.0.) THEN WRITE(lunout,*) 'this is shit' WRITE(lunout,*) varname,'(',i,j,k,') =',varmth(i,j,k) ENDIF var(i,j,k,imth)=varmth(i,j,k) ENDDO ENDDO ENDDO ENDDO STATUS = NF_CLOSE(NCID) IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in closing file',status END SUBROUTINE get_aero_fromfile