Program resultatbcsnow C C This is a program witch read the 3 hours filese with bc in snow C and snow data. It calculates total burden in snow and burden in snow C north of 65 degrees C C We will also calculate snow concentrations to be compared with C observations. C C The file also includes a subroutine for reading the netcdf file. Implicit None character*80 filename character*80 path Integer, parameter :: & times= 8, & lats=64, & lons=128 Integer :: i, j, k Integer :: ix, iy, it real time_vals(times) real niv_vals(10) real lat_vals(lats) real lon_vals(lons) real snowlayr(lons, lats, 10, times) real bcffcsnw(lons, lats, 10, times) real bcbiosnw(lons, lats, 10, times) C To calculate burden real :: jordradius real :: radianer real sumffc, sumbio, sumsno, Area, SUMA real sumffcN65, sumbioN65, sumsnoN65 C To calculate concentrations real SUMsnowlayr(lons, lats) real SUMbcffcsnw(lons, lats) real SUMbcbiosnw(lons, lats) real SUMconc(lons,lats) real SUMsnowlayr1(lons, lats) real SUMbcffcsnw1(lons, lats) real SUMbcbiosnw1(lons, lats) real SUMconc1(lons,lats) real SUMsnowlayr2(lons, lats) real SUMbcffcsnw2(lons, lats) real SUMbcbiosnw2(lons, lats) real SUMconc2(lons,lats) real SUMsnowlayr3(lons, lats) real SUMbcffcsnw3(lons, lats) real SUMbcbiosnw3(lons, lats) real SUMconc3(lons,lats) cMTL real SUMsnowlayr4(lons, lats) real SUMbcffcsnw4(lons, lats) real SUMbcbiosnw4(lons, lats) real SUMconc4(lons,lats) real SUMsnowlayr5(lons, lats) real SUMbcffcsnw5(lons, lats) real SUMbcbiosnw5(lons, lats) real SUMconc5(lons,lats) real SUMsnowlayr6(lons, lats) real SUMbcffcsnw6(lons, lats) real SUMbcbiosnw6(lons, lats) real SUMconc6(lons,lats) cMTL real SUMsnowlayrYEAR(lons, lats) real SUMbcffcsnwYEAR(lons, lats) real SUMbcbiosnwYEAR(lons, lats) real SUMconcYEAR(lons,lats) integer tellerxy(lons,lats) integer tellerxy1(lons,lats) integer tellerxy2(lons,lats) integer tellerxy3(lons,lats) integer tellerxyYEAR(lons,lats) logical fortsett integer tx, l cMTL integer tellerxy4(lons,lats) integer tellerxy5(lons,lats) integer tellerxy6(lons,lats) cMTL C For concentration at stations C Antall stasjoner Integer, parameter :: c & stant = 23+1 & stant = 21 character*3 sttext integer stasjnr(stant) ! stasjonsnummer integer staslon(stant) ! stasjonens lengdegrad integer staslat(stant) ! stasjonens breddegrad character*8 stasjname !stasjonsnavn integer rstat integer year, yr character*4 :: yeartext integer yeararray(2) integer day character*3 :: daytext character*80 fileout1 character*80 fileout2 character*80 fileout3 character*80 fileout4 character*20 project character*20 scen real :: temp real ::temp2 real ::temp3 real ::temp4 real ::bctemp real ::bcconc real ::snow C Input: choose year data yeararray/2008,2009/ scen = 'eclipse2008_100ice' project = '2013' DO yr = 1,1 year = 2008 !yeararray(yr) write(yeartext,FMT='(i4.4)') year path = '/div/enso/d10//blackcarbon2/PhD/' // & TRIM(project) // '/' // TRIM(scen) //'/snw/' C Output: C Burden: fileout1 = 'BCSNOW/' & // trim(scen) // & '/global_bcsnow_burden_' //yeartext // '.txt' open(21, FILE=fileout1) C Burden north of 65 degrees fileout2 = 'BCSNOW/' & // trim(scen) // & '/N65_bcsnow_burden_' //yeartext // '.txt' open(22, FILE=fileout2) C Mean concentration MAM fileout3 = 'BCSNOW/' & // trim(scen) // & '/MAM_bcconc_'//yeartext // '.txt' open(23, FILE=fileout3) C Mean concentration annual fileout4 = 'BCSNOW/' & // trim(scen) // & '/ANNUAL_bcconc_'//yeartext // '.txt' open(24, FILE=fileout4) C Mean concentration MAM fileout3 = 'BCSNOW/' & // trim(scen) // & '/MAR_bcconc_'//yeartext // '.txt' open(25, FILE=fileout3) C Mean concentration MAM fileout3 = 'BCSNOW/' & // trim(scen) // & '/APR_bcconc_'//yeartext // '.txt' open(26, FILE=fileout3) C Mean concentration MAM fileout3 = 'BCSNOW/' & // trim(scen) // & '/MAY_bcconc_'//yeartext // '.txt' open(27, FILE=fileout3) cMTL Added JJA, SON and DJF output C Mean concentration JJA fileout3 = 'BCSNOW/' & // trim(scen) // & '/JJA_bcconc_'//yeartext // '.txt' open(28, FILE=fileout3) C Mean concentration SON fileout3 = 'BCSNOW/' & // trim(scen) // & '/SON_bcconc_'//yeartext // '.txt' open(29, FILE=fileout3) C Mean concentration DJF fileout3 = 'BCSNOW/' & // trim(scen) // & '/DJF_bcconc_'//yeartext // '.txt' open(30, FILE=fileout3) cMTL C Station output C Open file with station names, lon and lat. open(20, FILE='stasjBCsnow.dat',IOSTAT = rstat) print*, rstat C Loop through number of stations DO i = 1, stant read(UNIT=20,FMT='(A8,I4, I4 )',IOSTAT = rstat) & stasjname, staslon(i),staslat(i) if(rstat.NE.0) then print*, rstat, 'STOP HEI' stop endif print*, stasjname , 'Hei' , & staslon(i), 'Hei',staslat(i) stasjnr(i) = i write(sttext,FMT='(i3.3)') stasjnr(i) fileout3 = 'BCSNOW/' & // trim(scen) // & '/BCsnow_stasj_' // sttext // & '_' // yeartext // '.txt' open(30+stasjnr(i), FILE=fileout3) write(30+stasjnr(i),*) stasjname ENDDO close(20) C**************************************** C To calculate burden, hight.... C Jordradius: (m) jordradius = 6371*1000 radianer = 1 * 3.14 /180 C**************************************** Print*, 'Starting!' ! teller = 0 tellerxy = 0 SUMsnowlayr = 0.0 SUMbcffcsnw = 0.0 SUMbcbiosnw = 0.0 SUMconc = 0.0 tellerxy1 = 0 SUMsnowlayr1 = 0.0 SUMbcffcsnw1 = 0.0 SUMbcbiosnw1 = 0.0 SUMconc1 = 0.0 tellerxy2 = 0 SUMsnowlayr2 = 0.0 SUMbcffcsnw2 = 0.0 SUMbcbiosnw2 = 0.0 SUMconc2 = 0.0 tellerxy3 = 0 SUMsnowlayr3 = 0.0 SUMbcffcsnw3 = 0.0 SUMbcbiosnw3 = 0.0 SUMconc3 = 0.0 ! conc = 0.0 tellerxyYEAR = 0 SUMsnowlayrYEAR = 0.0 SUMbcffcsnwYEAR = 0.0 SUMbcbiosnwYEAR = 0.0 SUMconcYEAR = 0.0 C Loop through days DO day = 1, 365 Write(daytext,FMT='(i3.3)') day C**************************************** C Read data from file: C Fossil Fuel filename = & TRIM(path) // 'bcsnow_'//yeartext//'_' & // daytext // '.nc' Print *, filename call readbc_from(filename, snowlayr, bcffcsnw, bcbiosnw, & time_vals, & niv_vals, lat_vals, lon_vals) C DO it = 1, times SUMA = 0.0 Area = 0 sumsno = 0 sumffc = 0 sumbio = 0 sumsnoN65 = 0 sumffcN65 = 0 sumbioN65 = 0 do ix = 1, lats Area = jordradius*jordradius* & cos(lat_vals(ix)*radianer) & *2.8*radianer*2.8*radianer c Lengdegrader do iy = 1, lons !Summing the area of the globe SUMA = SUMA + Area !Summing the mass of bc in each gridbox !kg/m2*m2 -> kg sumffc = sumffc + sum(bcffcsnw(iy,ix,:,it))*Area sumbio = sumbio + sum(bcbiosnw(iy,ix,:,it))*Area sumsno = sumsno + sum(snowlayr(iy,ix,:,it))*Area !Summing the mass of bc north of 65N IF(lat_vals(ix).GT.65) THEN sumffcN65 = sumffcN65 + & sum(bcffcsnw(iy,ix,:,it))*Area sumbioN65 = sumbioN65 + & sum(bcbiosnw(iy,ix,:,it))*Area sumsnoN65 = sumsnoN65 + & sum(snowlayr(iy,ix,:,it))*Area ENDIF ENDDO !lons ENDDO !lats Print*, 'Check area (m2): ' , SUMA Print*, 'Check 4piR^2: ' , 4*3.14* jordradius*jordradius ! Print*, 'SUM (g) ', SUM ! Print*, 'SUMN65 (g) ', SUMN65 !Write burden in kg to file Write(21,FMT='(3I7,2F14.2,F25.5)') & day, it , year,sumffc, sumbio, sumsno Write(22,FMT='(3I7,2F14.2,F25.5)') & day, it , year,sumffcN65, sumbioN65, sumsnoN65 C Stasjonsoutput: DO i = 1, stant print*, stasjnr(i), & staslon(i),staslat(i) Write(30+stasjnr(i),FMT='(3I6,30E14.7)') day, it , & year, & snowlayr(staslon(i),staslat(i),:,it), & bcbiosnw(staslon(i),staslat(i),:,it), & bcffcsnw(staslon(i),staslat(i),:,it) IF(day.eq.1) then Print*, 'Station ', stasjnr(i) Print*, 'Location Lat' , lat_vals(staslat(i)) Print*, 'Location Lon' , lon_vals(staslon(i)) ENDIF ENDDO ENDDO !times ! MEAN CONCENTRATION CALCULATION PERIODE IF (day.ge.60.and.day.le.151) THEN do it = 1, times do ix = 1, lats do iy = 1, lons temp = 0.0 temp2 = 0.0 temp3 = 0.0 temp4 = 0.0 bctemp = 0.0 bcconc = 0.0 snow = 0.0 !IF hole snowlayer some thickness if(sum(snowlayr(iy,ix,:,it)).GE.0.001) then fortsett = .true. tx = 10 do while (fortsett) if(sum(snowlayr(iy,ix,1:tx,it)).gt.0.005)then tx = tx - 1 fortsett = .true. if(tx.eq.1) then fortsett = .false. C Summen opp til laget < 0.005 temp = sum(snowlayr(iy,ix,1:tx,it)) C Størrelsen på laget over temp2 = snowlayr(iy,ix,tx+1,it) C Snømengden vi mangler for å få 0.005 temp3 = 0.005 - temp C Andel vi mangler snø: temp4 = temp3/temp2 C Samme andel soot bctemp = temp4 & *(bcffcsnw(iy,ix,tx+1,it)+bcbiosnw(iy,ix,tx+1,it)) endif else ! Less than 0.005 fortsett = .false. temp = sum(snowlayr(iy,ix,1:tx,it)) if(tx == 10) then bctemp = 0.0 temp3 = 0.0 else temp2 = snowlayr(iy,ix,tx+1,it) temp3 = 0.005 - temp temp4 = temp3/temp2 bctemp = temp4 & *(bcffcsnw(iy,ix,tx+1,it)+bcbiosnw(iy,ix,tx+1,it)) endif endif enddo snow = sum(snowlayr(iy,ix,1:tx,it))+temp3 bcconc = (sum(bcffcsnw(iy,ix,1:tx,it))+ & sum(bcbiosnw(iy,ix,1:tx,it))+bctemp)/snow SUMsnowlayr(iy,ix) = SUMsnowlayr(iy,ix) + & snow C SUMbcffcsnw(iy,ix) = SUMbcffcsnw(iy,ix) + C & sum(bcffcsnw(iy,ix,1:tx,it)) C C C SUMbcbiosnw(iy,ix) = SUMbcbiosnw(iy,ix) + C & sum(bcbiosnw(iy,ix,1:tx,it)) C C SUMconc(iy,ix) = SUMconc(iy,ix) + & bcconc tellerxy(iy,ix) = tellerxy(iy,ix) + 1 endif enddo enddo enddo ! teller = teller + 1 ENDIF ! MEAN CONCENTRATION CALCULATION MAR IF (day.ge.60.and.day.le.90) THEN do it = 1, times do ix = 1, lats do iy = 1, lons temp = 0.0 temp2 = 0.0 temp3 = 0.0 temp4 = 0.0 bctemp = 0.0 bcconc = 0.0 snow = 0.0 !IF hole snowlayer some thickness if(sum(snowlayr(iy,ix,:,it)).GE.0.001) then fortsett = .true. tx = 10 do while (fortsett) if(sum(snowlayr(iy,ix,1:tx,it)).gt.0.005)then tx = tx - 1 fortsett = .true. if(tx.eq.1) then fortsett = .false. C Summen opp til laget < 0.005 temp = sum(snowlayr(iy,ix,1:tx,it)) C Størrelsen på laget over temp2 = snowlayr(iy,ix,tx+1,it) C Snømengden vi mangler for å få 0.005 temp3 = 0.005 - temp C Andel vi mangler snø: temp4 = temp3/temp2 C Samme andel soot bctemp = temp4 & *(bcffcsnw(iy,ix,tx+1,it)+bcbiosnw(iy,ix,tx+1,it)) endif else ! Less than 0.005 fortsett = .false. temp = sum(snowlayr(iy,ix,1:tx,it)) if(tx == 10) then bctemp = 0.0 temp3 = 0.0 else temp2 = snowlayr(iy,ix,tx+1,it) temp3 = 0.005 - temp temp4 = temp3/temp2 bctemp = temp4 & *(bcffcsnw(iy,ix,tx+1,it)+bcbiosnw(iy,ix,tx+1,it)) endif endif enddo snow = sum(snowlayr(iy,ix,1:tx,it))+temp3 bcconc = (sum(bcffcsnw(iy,ix,1:tx,it))+ & sum(bcbiosnw(iy,ix,1:tx,it))+bctemp)/snow SUMsnowlayr1(iy,ix) = SUMsnowlayr1(iy,ix) + & snow C SUMbcffcsnw1(iy,ix) = SUMbcffcsnw1(iy,ix) + C & sum(bcffcsnw(iy,ix,1:tx,it)) C C C SUMbcbiosnw1(iy,ix) = SUMbcbiosnw1(iy,ix) + C & sum(bcbiosnw(iy,ix,1:tx,it)) C C SUMconc1(iy,ix) = SUMconc1(iy,ix) + & bcconc tellerxy1(iy,ix) = tellerxy1(iy,ix) + 1 endif enddo enddo enddo ! teller = teller + 1 ENDIF ! MEAN CONCENTRATION CALCULATION APRIL IF (day.ge.91.and.day.le.120) THEN do it = 1, times do ix = 1, lats do iy = 1, lons temp = 0.0 temp2 = 0.0 temp3 = 0.0 temp4 = 0.0 bctemp = 0.0 bcconc = 0.0 snow = 0.0 !IF hole snowlayer some thickness if(sum(snowlayr(iy,ix,:,it)).GE.0.001) then fortsett = .true. tx = 10 do while (fortsett) if(sum(snowlayr(iy,ix,1:tx,it)).gt.0.005)then tx = tx - 1 fortsett = .true. if(tx.eq.1) then fortsett = .false. C Summen opp til laget < 0.005 temp = sum(snowlayr(iy,ix,1:tx,it)) C Størrelsen på laget over temp2 = snowlayr(iy,ix,tx+1,it) C Snømengden vi mangler for å få 0.005 temp3 = 0.005 - temp C Andel vi mangler snø: temp4 = temp3/temp2 C Samme andel soot bctemp = temp4 & *(bcffcsnw(iy,ix,tx+1,it)+bcbiosnw(iy,ix,tx+1,it)) endif else ! Less than 0.005 fortsett = .false. temp = sum(snowlayr(iy,ix,1:tx,it)) if(tx == 10) then bctemp = 0.0 temp3 = 0.0 else temp2 = snowlayr(iy,ix,tx+1,it) temp3 = 0.005 - temp temp4 = temp3/temp2 bctemp = temp4 & *(bcffcsnw(iy,ix,tx+1,it)+bcbiosnw(iy,ix,tx+1,it)) endif endif enddo snow = sum(snowlayr(iy,ix,1:tx,it))+temp3 bcconc = (sum(bcffcsnw(iy,ix,1:tx,it))+ & sum(bcbiosnw(iy,ix,1:tx,it))+bctemp)/snow SUMsnowlayr2(iy,ix) = SUMsnowlayr2(iy,ix) + & snow C SUMbcffcsnw2(iy,ix) = SUMbcffcsnw2(iy,ix) + C & sum(bcffcsnw(iy,ix,1:tx,it)) C C C SUMbcbiosnw2(iy,ix) = SUMbcbiosnw2(iy,ix) + C & sum(bcbiosnw(iy,ix,1:tx,it)) C C SUMconc2(iy,ix) = SUMconc2(iy,ix) + & bcconc tellerxy2(iy,ix) = tellerxy2(iy,ix) + 1 endif enddo enddo enddo ! teller = teller + 1 ENDIF ! MEAN CONCENTRATION CALCULATION MAY IF (day.ge.121.and.day.le.151) THEN do it = 1, times do ix = 1, lats do iy = 1, lons temp = 0.0 temp2 = 0.0 temp3 = 0.0 temp4 = 0.0 bctemp = 0.0 bcconc = 0.0 snow = 0.0 !IF hole snowlayer some thickness if(sum(snowlayr(iy,ix,:,it)).GE.0.001) then fortsett = .true. tx = 10 do while (fortsett) if(sum(snowlayr(iy,ix,1:tx,it)).gt.0.005)then tx = tx - 1 fortsett = .true. if(tx.eq.1) then fortsett = .false. C Summen opp til laget < 0.005 temp = sum(snowlayr(iy,ix,1:tx,it)) C Størrelsen på laget over temp2 = snowlayr(iy,ix,tx+1,it) C Snømengden vi mangler for å få 0.005 temp3 = 0.005 - temp C Andel vi mangler snø: temp4 = temp3/temp2 C Samme andel soot bctemp = temp4 & *(bcffcsnw(iy,ix,tx+1,it)+bcbiosnw(iy,ix,tx+1,it)) endif else ! Less than 0.005 fortsett = .false. temp = sum(snowlayr(iy,ix,1:tx,it)) if(tx == 10) then bctemp = 0.0 temp3 = 0.0 else temp2 = snowlayr(iy,ix,tx+1,it) temp3 = 0.005 - temp temp4 = temp3/temp2 bctemp = temp4 & *(bcffcsnw(iy,ix,tx+1,it)+bcbiosnw(iy,ix,tx+1,it)) endif endif enddo snow = sum(snowlayr(iy,ix,1:tx,it))+temp3 bcconc = (sum(bcffcsnw(iy,ix,1:tx,it))+ & sum(bcbiosnw(iy,ix,1:tx,it))+bctemp)/snow SUMsnowlayr3(iy,ix) = SUMsnowlayr3(iy,ix) + & snow C SUMbcffcsnw3(iy,ix) = SUMbcffcsnw1(iy,ix) + C & sum(bcffcsnw(iy,ix,1:tx,it)) C C C SUMbcbiosnw3(iy,ix) = SUMbcbiosnw1(iy,ix) + C & sum(bcbiosnw(iy,ix,1:tx,it)) C C SUMconc3(iy,ix) = SUMconc3(iy,ix) + & bcconc tellerxy3(iy,ix) = tellerxy3(iy,ix) + 1 endif enddo enddo enddo ! teller = teller + 1 ENDIF cMTL ! MEAN CONCENTRATION CALCULATION JJA IF (day.ge.152.and.day.le.243) THEN do it = 1, times do ix = 1, lats do iy = 1, lons temp = 0.0 temp2 = 0.0 temp3 = 0.0 temp4 = 0.0 bctemp = 0.0 bcconc = 0.0 snow = 0.0 !IF hole snowlayer some thickness if(sum(snowlayr(iy,ix,:,it)).GE.0.001) then fortsett = .true. tx = 10 do while (fortsett) if(sum(snowlayr(iy,ix,1:tx,it)).gt.0.005)then tx = tx - 1 fortsett = .true. if(tx.eq.1) then fortsett = .false. C Summen opp til laget < 0.005 temp = sum(snowlayr(iy,ix,1:tx,it)) C Størrelsen på laget over temp2 = snowlayr(iy,ix,tx+1,it) C Snømengden vi mangler for å få 0.005 temp3 = 0.005 - temp C Andel vi mangler snø: temp4 = temp3/temp2 C Samme andel soot bctemp = temp4 & *(bcffcsnw(iy,ix,tx+1,it)+bcbiosnw(iy,ix,tx+1,it)) endif else ! Less than 0.005 fortsett = .false. temp = sum(snowlayr(iy,ix,1:tx,it)) if(tx == 10) then bctemp = 0.0 temp3 = 0.0 else temp2 = snowlayr(iy,ix,tx+1,it) temp3 = 0.005 - temp temp4 = temp3/temp2 bctemp = temp4 & *(bcffcsnw(iy,ix,tx+1,it)+bcbiosnw(iy,ix,tx+1,it)) endif endif enddo snow = sum(snowlayr(iy,ix,1:tx,it))+temp3 bcconc = (sum(bcffcsnw(iy,ix,1:tx,it))+ & sum(bcbiosnw(iy,ix,1:tx,it))+bctemp)/snow SUMsnowlayr4(iy,ix) = SUMsnowlayr4(iy,ix) + & snow C SUMbcffcsnw4(iy,ix) = SUMbcffcsnw4(iy,ix) + C & sum(bcffcsnw(iy,ix,1:tx,it)) C C C SUMbcbiosnw4(iy,ix) = SUMbcbiosnw4(iy,ix) + C & sum(bcbiosnw(iy,ix,1:tx,it)) C C SUMconc4(iy,ix) = SUMconc4(iy,ix) + & bcconc tellerxy4(iy,ix) = tellerxy4(iy,ix) + 1 endif enddo enddo enddo ! teller = teller + 1 ENDIF ! MEAN CONCENTRATION CALCULATION SON IF (day.ge.244.and.day.le.334) THEN do it = 1, times do ix = 1, lats do iy = 1, lons temp = 0.0 temp2 = 0.0 temp3 = 0.0 temp4 = 0.0 bctemp = 0.0 bcconc = 0.0 snow = 0.0 !IF hole snowlayer some thickness if(sum(snowlayr(iy,ix,:,it)).GE.0.001) then fortsett = .true. tx = 10 do while (fortsett) if(sum(snowlayr(iy,ix,1:tx,it)).gt.0.005)then tx = tx - 1 fortsett = .true. if(tx.eq.1) then fortsett = .false. C Summen opp til laget < 0.005 temp = sum(snowlayr(iy,ix,1:tx,it)) C Størrelsen på laget over temp2 = snowlayr(iy,ix,tx+1,it) C Snømengden vi mangler for å få 0.005 temp3 = 0.005 - temp C Andel vi mangler snø: temp4 = temp3/temp2 C Samme andel soot bctemp = temp4 & *(bcffcsnw(iy,ix,tx+1,it)+bcbiosnw(iy,ix,tx+1,it)) endif else ! Less than 0.005 fortsett = .false. temp = sum(snowlayr(iy,ix,1:tx,it)) if(tx == 10) then bctemp = 0.0 temp3 = 0.0 else temp2 = snowlayr(iy,ix,tx+1,it) temp3 = 0.005 - temp temp4 = temp3/temp2 bctemp = temp4 & *(bcffcsnw(iy,ix,tx+1,it)+bcbiosnw(iy,ix,tx+1,it)) endif endif enddo snow = sum(snowlayr(iy,ix,1:tx,it))+temp3 bcconc = (sum(bcffcsnw(iy,ix,1:tx,it))+ & sum(bcbiosnw(iy,ix,1:tx,it))+bctemp)/snow SUMsnowlayr5(iy,ix) = SUMsnowlayr5(iy,ix) + & snow C SUMbcffcsnw5(iy,ix) = SUMbcffcsnw5(iy,ix) + C & sum(bcffcsnw(iy,ix,1:tx,it)) C C C SUMbcbiosnw5(iy,ix) = SUMbcbiosnw5(iy,ix) + C & sum(bcbiosnw(iy,ix,1:tx,it)) C C SUMconc5(iy,ix) = SUMconc5(iy,ix) + & bcconc tellerxy5(iy,ix) = tellerxy5(iy,ix) + 1 endif enddo enddo enddo ! teller = teller + 1 ENDIF ! MEAN CONCENTRATION CALCULATION DJF IF (day.ge.1.and.day.le.59 .or. & day .ge. 335 .and. day .le. 365) THEN do it = 1, times do ix = 1, lats do iy = 1, lons temp = 0.0 temp2 = 0.0 temp3 = 0.0 temp4 = 0.0 bctemp = 0.0 bcconc = 0.0 snow = 0.0 !IF hole snowlayer some thickness if(sum(snowlayr(iy,ix,:,it)).GE.0.001) then fortsett = .true. tx = 10 do while (fortsett) if(sum(snowlayr(iy,ix,1:tx,it)).gt.0.005)then tx = tx - 1 fortsett = .true. if(tx.eq.1) then fortsett = .false. C Summen opp til laget < 0.005 temp = sum(snowlayr(iy,ix,1:tx,it)) C Størrelsen på laget over temp2 = snowlayr(iy,ix,tx+1,it) C Snømengden vi mangler for å få 0.005 temp3 = 0.005 - temp C Andel vi mangler snø: temp4 = temp3/temp2 C Samme andel soot bctemp = temp4 & *(bcffcsnw(iy,ix,tx+1,it)+bcbiosnw(iy,ix,tx+1,it)) endif else ! Less than 0.005 fortsett = .false. temp = sum(snowlayr(iy,ix,1:tx,it)) if(tx == 10) then bctemp = 0.0 temp3 = 0.0 else temp2 = snowlayr(iy,ix,tx+1,it) temp3 = 0.005 - temp temp4 = temp3/temp2 bctemp = temp4 & *(bcffcsnw(iy,ix,tx+1,it)+bcbiosnw(iy,ix,tx+1,it)) endif endif enddo snow = sum(snowlayr(iy,ix,1:tx,it))+temp3 bcconc = (sum(bcffcsnw(iy,ix,1:tx,it))+ & sum(bcbiosnw(iy,ix,1:tx,it))+bctemp)/snow SUMsnowlayr6(iy,ix) = SUMsnowlayr6(iy,ix) + & snow C SUMbcffcsnw6(iy,ix) = SUMbcffcsnw6(iy,ix) + C & sum(bcffcsnw(iy,ix,1:tx,it)) C C C SUMbcbiosnw6(iy,ix) = SUMbcbiosnw6(iy,ix) + C & sum(bcbiosnw(iy,ix,1:tx,it)) C C SUMconc6(iy,ix) = SUMconc6(iy,ix) + & bcconc tellerxy6(iy,ix) = tellerxy6(iy,ix) + 1 endif enddo enddo enddo ! teller = teller + 1 ENDIF cMTL ! MEAN CONCENTRATION CALCULATION WHOLE YEAR do it = 1, times do ix = 1, lats do iy = 1, lons temp = 0.0 temp2 = 0.0 temp3 = 0.0 temp4 = 0.0 bctemp = 0.0 bcconc = 0.0 snow = 0.0 !IF hole snowlayer some thickness if(sum(snowlayr(iy,ix,:,it)).GE.0.001) then fortsett = .true. tx = 10 do while (fortsett) if(sum(snowlayr(iy,ix,1:tx,it)).gt.0.005)then tx = tx - 1 fortsett = .true. if(tx.eq.1) then fortsett = .false. C Summen opp til laget < 0.005 temp = sum(snowlayr(iy,ix,1:tx,it)) C Størrelsen på laget over temp2 = snowlayr(iy,ix,tx+1,it) C Snømengden vi mangler for å få 0.005 temp3 = 0.005 - temp C Andel vi mangler snø: temp4 = temp3/temp2 C Samme andel soot bctemp = temp4 & *(bcffcsnw(iy,ix,tx+1,it)+bcbiosnw(iy,ix,tx+1,it)) endif else ! Less than 0.005 fortsett = .false. temp = sum(snowlayr(iy,ix,1:tx,it)) if(tx == 10) then bctemp = 0.0 temp3 = 0.0 else temp2 = snowlayr(iy,ix,tx+1,it) temp3 = 0.005 - temp temp4 = temp3/temp2 bctemp = temp4 & *(bcffcsnw(iy,ix,tx+1,it)+bcbiosnw(iy,ix,tx+1,it)) endif endif enddo snow = sum(snowlayr(iy,ix,1:tx,it))+temp3 bcconc = (sum(bcffcsnw(iy,ix,1:tx,it))+ & sum(bcbiosnw(iy,ix,1:tx,it))+bctemp)/snow SUMsnowlayrYEAR(iy,ix) = & SUMsnowlayrYEAR(iy,ix) + & snow C SUMbcffcsnw1(iy,ix) = SUMbcffcsnw1(iy,ix) + C & sum(bcffcsnw(iy,ix,1:tx,it)) C C C SUMbcbiosnw1(iy,ix) = SUMbcbiosnw1(iy,ix) + C & sum(bcbiosnw(iy,ix,1:tx,it)) C C SUMconcYEAR(iy,ix) = SUMconcYEAR(iy,ix) + & bcconc tellerxyYEAR(iy,ix) = tellerxyYEAR(iy,ix) + 1 endif enddo enddo enddo ! teller = teller + 1 ENDDO !day do k = 1,lats do l = 1,lons write(23,*) l,k, SUMconc(l,k)/tellerxy(l,k), & SUMsnowlayr(l,k)/tellerxy(l,k) enddo enddo do k = 1,lats do l = 1,lons write(24,*) l,k, SUMconcYEAR(l,k)/tellerxyYEAR(l,k), & SUMsnowlayrYEAR(l,k)/tellerxyYEAR(l,k) enddo enddo do k = 1,lats do l = 1,lons write(25,*) l,k, SUMconc1(l,k)/tellerxy1(l,k), & SUMsnowlayr1(l,k)/tellerxy1(l,k) enddo enddo do k = 1,lats do l = 1,lons write(26,*) l,k, SUMconc2(l,k)/tellerxy2(l,k), & SUMsnowlayr2(l,k)/tellerxy2(l,k) enddo enddo do k = 1,lats do l = 1,lons write(27,*) l,k, SUMconc3(l,k)/tellerxy3(l,k), & SUMsnowlayr3(l,k)/tellerxy3(l,k) enddo enddo cMTL do k = 1,lats do l = 1,lons write(28,*) l,k, SUMconc4(l,k)/tellerxy4(l,k), & SUMsnowlayr4(l,k)/tellerxy4(l,k) enddo enddo do k = 1,lats do l = 1,lons write(29,*) l,k, SUMconc5(l,k)/tellerxy5(l,k), & SUMsnowlayr5(l,k)/tellerxy5(l,k) enddo enddo do k = 1,lats do l = 1,lons write(30,*) l,k, SUMconc6(l,k)/tellerxy6(l,k), & SUMsnowlayr6(l,k)/tellerxy6(l,k) enddo enddo cMTL close(21) close(22) close(23) close(24) close(25) close(26) close(27) cMTL close(28) close(29) close(30) cMTL DO i = 1, stant close(30+stasjnr(i)) ENDDO ENDDO ! YEAR END C Reading bc and snow data from netcdf-file Subroutine readbc_from(filename, & snowlayr, bcffcsnw, bcbiosnw, time_vals, & niv_vals, lat_vals, lon_vals) use netcdf Implicit None character*80 filename integer err integer fileid integer dimid_lon,dimid_lat,dimid_time,dimid_niv integer varid_lon, varid_lat, varid_niv, varid_time integer varid_snowlayr, varid_bcffcsnw, varid_bcbiosnw Integer, parameter :: & times= 8, & lats=64, & lons=128 real time_vals(times) real niv_vals(10) real lat_vals(lats) real lon_vals(lons) integer start_1D(1),count_1D(1) integer start_4D(4),count_4D(4) real snowlayr(lons, lats, 10, times) real bcffcsnw(lons, lats, 10, times) real bcbiosnw(lons, lats, 10, times) integer t,k,i,j C Open the netCDF file for reading err=nf90_open( trim(filename), nf90_nowrite, fileid) Print*,'Opened the file!' , fileid c$$$C Get dimension IDs Print*, 'Get dimensionIDs' err=NF90_INQ_DIMID (fileid,'time',dimid_time) c print*,'time',dimid_time err=NF90_INQ_DIMID (fileid,'niv',dimid_niv) c print*,'lev',dimid_lev err=NF90_INQ_DIMID (fileid,'lat',dimid_lat) c print*,'lat',dimid_lat err=NF90_INQ_DIMID (fileid,'lon',dimid_lon) c print*,'lon', dimid_lon c$$$C Get variable IDs Print *, 'Get variable IDs' err=NF90_INQ_VARID (fileid,'time',varid_time) c print*,'varid_time',varid_time err=NF90_INQ_VARID (fileid,'niv',varid_niv) c print*,'varied_niv',varid_niv err=NF90_INQ_VARID (fileid,'lat',varid_lat) c print*,'varid_lat',varid_lat err=NF90_INQ_VARID (fileid,'lon',varid_lon) c print*,'varid_lon',varid_lon err=NF90_INQ_VARID (fileid,'snowlayr',varid_snowlayr) print*,'varid_snowlayr',varid_snowlayr err=NF90_INQ_VARID (fileid,'bcffcsnw',varid_bcffcsnw) print*,'varid_bcffcsnw',varid_bcffcsnw err=NF90_INQ_VARID (fileid,'bcbiosnw',varid_bcbiosnw) print*,'varid_bcbiosnw',varid_bcbiosnw C Get contents of dimension variables start_1D=1 count_1D=lons err=NF90_GET_VAR(fileid,varid_lon,lon_vals,start_1D,count_1D) c Print*, 'Err: ', err start_1D=1 count_1D=lats err=NF90_GET_VAR(fileid,varid_lat,lat_vals,start_1D,count_1D) c Print*, 'Err: ', err start_1D=1 count_1D=10 err=NF90_GET_VAR(fileid,varid_niv,niv_vals,start_1D,count_1D) c Print*, 'Err: ', err start_1D=1 count_1D=times err=NF90_GET_VAR(fileid,varid_time,time_vals,start_1D,count_1D) c Print*, 'Err: ', err Print *, 'Get the main data' start_4D(1)=1 start_4D(2)=1 start_4D(3)=1 start_4D(4)=1 ! outermost (ncdump) index first ! count_4D(1)=lons count_4D(2)=lats count_4D(3)=10 count_4D(4)=times ! innermost (ncdump) index last ! c Print*, start_4D c Print*, count_4D c Print*, varid_BCPHOB err=NF90_GET_VAR(fileid,varid_snowlayr,snowlayr,start_4D,count_4D) print*,'err for PS:',err start_4D(1)=1 start_4D(2)=1 start_4D(3)=1 start_4D(4)=1 ! outermost (ncdump) index first ! count_4D(1)=lons count_4D(2)=lats count_4D(3)=10 count_4D(4)=times ! innermost (ncdump) index last ! c Print*, count_4D c Print*, varid_BCPHIL err=NF90_GET_VAR(fileid,varid_bcffcsnw,bcffcsnw,start_4D,count_4D) print*,'err for PS:',err start_4D(1)=1 start_4D(2)=1 start_4D(3)=1 start_4D(4)=1 ! outermost (ncdump) index first ! count_4D(1)=lons count_4D(2)=lats count_4D(3)=10 count_4D(4)=times ! innermost (ncdump) index last ! c Print*, count_4D c Print*, varid_BCPHIL err=NF90_GET_VAR(fileid,varid_bcbiosnw,bcbiosnw,start_4D,count_4D) print*,'err for PS:',err err=NF90_CLOSE(fileid) Print *, 'Closed file!', err end subroutine readbc_from