c to do: c take init_tr out of comm_tracer.h, put inittr(MAXTRAC) in comm_new.h c allocate ntrac_tot of p_trclim in model_memory, if needed, otherwise, c delete from array_init #define c_str(s) ('s\0') c------------------------------------------------------------ subroutine hfx_pert_init c------------------------------------------------------------ include 'comm_para.h' include 'comm_new.h' include 'comm_data.h' include 'comm_tracer.h' real inp_flt, inp_days logical inp_def dimension flt(100) ihfprt = inp_int(c_str(Hflx_prt),0) if (ihfprt .gt. 0) then hfprt_amp = inp_flt(c_str(Hflx_prt_amp),15.0) hfprt_lat = inp_flt(c_str(Hflx_prt_lat),20.0) endif return end c------------------------------------------------------------ subroutine wind_stir_init c------------------------------------------------------------ include 'comm_para.h' include 'comm_new.h' include 'comm_data.h' include 'comm_tracer.h' real inp_flt, inp_days logical inp_def dimension flt(100) windm0 = inp_flt(c_str(m0_WndStir),0.6e-3) windz0 = inp_flt(c_str(z0_WndStir),50.0) return end c------------------------------------------------------------ subroutine tracer_input(npt,nz,ntimes,nstart,nstep) c------------------------------------------------------------ implicit real(a-h,o-z),integer(i-n) include 'comm_para.h' include 'comm_new.h' include 'comm_data.h' include 'comm_tracer.h' real inp_flt, inp_days logical inp_def dimension flt(100),xga(MAXTR),yga(MAXTR),dga(MAXTR),rga(MAXTR) common /gauss_trac/ xga, yga, dga, rga call array_init(npt,nz) delt = inp_days (c_str(Time_step), 1./24.) stpd = 1./delt ntrcont = int(stpd*inp_days(c_str(Tr_int_step),15.)) init_tr = inp_int(c_str(Tr_init),0) if (init_tr.ne.0) then print*,'+Tr_init is now obsolete' print*,'use the +Trac_init [array], instead' stop endif ifilt_tr = inp_int(c_str(Tr_filt),1) iforc_tr = inp_int(c_str(Tr_forcing),0) if (iforc_tr.ne.0) then print*,'+TR_forcing is now obsolete' print*,'use the +Trac_forc [array], instead' stop endif nret = inp_iarr (c_str(Trac_forc), ntrac_tot, forctr, forctr) do itrac = 1, ntrac if (forctr(itrac).eq.5 ) forctr(itrac+1) = 0 if (forctr(itrac).eq.6.or.forctr(itrac).eq.8) then forctr(itrac+1) = 0 forctr(itrac+2) = 0 endif if (forctr(itrac).eq.9) then do it = 1, 9 forctr(itrac+it) = 0 enddo endif if (forctr(itrac).eq.7 ) forctr(itrac+1) = 0 if (forctr(itrac).eq.53) forctr(itrac+1) = 0 enddo icl_tr = inp_int(c_str(Clim_Tr),0) ipp_tr = inp_int(c_str(Tr_pp),0) itanom_init = inp_int(c_str(Tanom_init),0) itrac_cons = inp_int(c_str(Tr_cons_out),0) igas_ex = inp_int(c_str(Gas_Ex),0) c-*-*-*-*-*-* c-*-*-*-*-*-* c-*-*-*-*-*-* c-*-*-*-*-*-* c BEEF THIS SECTION UP !!!!!!! if (forctr(1).eq.6 .and. ntrac.ne.3) then write(*,*) " forcing option inconsistent with ntrac" endif c-*-*-*-*-*-* c-*-*-*-*-*-* c-*-*-*-*-*-* c-*-*-*-*-*-* c-*-*-*-*-*-* i = inp_rarr(c_str(Tr_gs_lon), ntrac, xga, xga) j = inp_rarr(c_str(Tr_gs_lat), ntrac, yga, yga) k = inp_rarr(c_str(Tr_gs_dp), ntrac, dga, dga) l = inp_rarr(c_str(Tr_gs_rad), ntrac, rga, rga) do jtr = 1, ntrac if (inittr(jtr).eq.21) then if (i.le.jtr.or.j.le.jtr.or.k.le.jtr.or.l.le.jtr) then print*,'Not enough data in +Tr_gs_* arrays' stop endif endif enddo c c for c14 with constant exchange, c read co2 gas exchange flux in units c of moles CO2/m2/yr from input file c if (forctr(1).ge.1 .or. forctr(1).le.10) then co2geflx = inp_flt(c_str(CO2_gflx),20.) factor_c14 = co2geflx/(100.*24.*365.*3600.) endif if (igas_ex.eq.1 .or. igas_ex.eq.2) then n_atm = inp_str(c_str(Tr_atm_forc),'none',fatf) call odb_open(idf_tr(1),fatf(1:n_atm),0) call odb_rddm(idf_tr(1),'T',nt_tratm) call odb_rddm(idf_tr(1),'LAT',nlat_tratm) call mem_alloc(p_tr_atm, nt_tratm*nlat_tratm,2,'tr_atm') call mem_alloc(p_tr_tgrid,nt_tratm,2,'tr_tgrid') call mem_alloc(p_tr_latgrid,nlat_tratm,2,'tr_latgrid') call odb_rdgr(idf_tr(1),'T',nt_tratm,tr_tgrid) call odb_rdgr(idf_tr(1),'LAT',nlat_tratm,tr_latgrid) if ((iforc_tr.ge.1).and.(iforc_tr.le.10)) then call odb_rdvar(idf_tr(1),'c14',tr_atm) endif endif return end c---------------------------------------------------- subroutine tritium_init(npt,nz,nstart,nstep) c---------------------------------------------------- include 'comm_new.h' include 'comm_data.h' include 'comm_tracer.h' character c70*70, c5*5 data ireawr /0/ data ireadon /0/ pi = asin(1.)*2. c c Read in doney's arrays c f1,f2,am,phz c n_rdonf1 = inp_str(c_str(Tr_donf1) ,'none',fdonf1) n_rdonf2 = inp_str(c_str(Tr_donf2) ,'none',fdonf2) n_rdonam = inp_str(c_str(Tr_donam) ,'none',fdonam) n_rdonphz = inp_str(c_str(Tr_donphz) ,'none',fdonphz) call odb_open(idf_donf1, fdonf1(1:n_rdonf1),0) call odb_open(idf_donf2, fdonf2(1:n_rdonf2),0) call odb_open(idf_donam, fdonam(1:n_rdonam),0) call odb_open(idf_donphz,fdonphz(1:n_rdonphz),0) call data_on_model_grid(idf_donf1, ldonf1, 'F1') call read_zt(idf_donf1,ldonf1,npt,1,1,'F1',donf1(1)) call data_on_model_grid(idf_donf2, ldonf2, 'F2') call read_zt(idf_donf2,ldonf2,npt,1,1,'F2',donf2(1)) call data_on_model_grid(idf_donam, ldonam, 'AM') call read_zt(idf_donam,ldonam,npt,1,1,'AM',donam(1)) call data_on_model_grid(idf_donf1, ldonphz, 'F1') call read_zt(idf_donphz,ldonphz,npt,1,1,'PHZ',donphz(1)) c c read in atmospheric data: c relative humidity(m/yr) c evaporation(m/yr) c precipitation(m/yr) c n_rdonrh = inp_str(c_str(Tr_donrh) ,'none',fdonrh) n_rdonevp = inp_str(c_str(Tr_donevp) ,'none',fdonevp) n_rdonprcp = inp_str(c_str(Tr_donprcp),'none',fdonprcp) n_rdonabwn = inp_str(c_str(Tr_donabwn),'none',fdonabwn) call odb_open(idf_donrh, fdonrh(1:n_rdonrh),0) call odb_open(idf_donevp, fdonevp(1:n_rdonevp),0) call odb_open(idf_donprcp,fdonprcp(1:n_rdonprcp),0) call odb_open(idf_donabwn,fdonabwn(1:n_rdonabwn),0) call odb_rddm(idf_donevp, 'T', nevap) call mem_alloc(p_tdoney,nevap,2,'tdoney') call odb_rdgr(idf_donevp,'T',nevap, tdoney) call it_catch(nevap,tdoney,nstart,it1,it2,tscl) idoney = it2 call data_on_model_grid(idf_donevp, levap, 'evap') call read_zt(idf_donevp,levap,npt,1,it1,'evap',evap(1)) call read_zt(idf_donevp,levap,npt,1,it2,'evap',evap(1+npt)) call data_on_model_grid(idf_donprcp, lprecip, 'precip') call read_zt(idf_donprcp,lprecip,npt,1,it1,'precip',precip(1)) call read_zt(idf_donprcp,lprecip,npt,1,it2,'precip',precip(1+npt)) call data_on_model_grid(idf_donrh, lrhum, 'relhum') call read_zt(idf_donrh,lrhum,npt,1,it1,'relhum',relhum(1)) call read_zt(idf_donrh,lrhum,npt,1,it2,'relhum',relhum(1+npt)) call data_on_model_grid(idf_donabwn, labwn, 'abswin') call read_zt(idf_donabwn,labwn,npt,1,it1,'abswin',abswin(1)) call read_zt(idf_donabwn,labwn,npt,1,it2,'abswin',abswin(1+npt)) do i=1,npt evap(i+2*npt)=evap(i)+tscl*(evap(i+npt)-evap(i)) precip(i+2*npt)=precip(i)+tscl*(precip(i+npt)-precip(i)) relhum(i+2*npt)=relhum(i)+tscl*(relhum(i+npt)-relhum(i)) abswin(i+2*npt)=abswin(i)+tscl*(abswin(i+npt)-abswin(i)) enddo c BEGIN STUFF FROM CHRISTOPH'S CODE: c (read ascii files) open(15,file='/home/keithr/MODEL/tritium_source_wr', * access='sequential', form='formatted') open(19,file='/home/keithr/MODEL/doney_input/factor_scores.dat', * form='formatted',status='unknown') c---compute monthly rate of tritium input according c to roether and weiss (1980), constant source over c one year; i.e. the source is compute newly once c a year c c read time curves c print*,' !!!!!!!!!! ireawr ',ireawr if(ireawr.eq.0)then rewind 15 read(15,'(a70)')c70 do l=1,newr read(15,'(f7.1,6(3x,f7.1))') * souryr(l),cp50n(l),cp50s(l),cr50n(l), * sp50n(l),sp50s(l),sr50n(l) c write(6,'(i2,1x,f7.1,6(3x,f7.1))') c * l,souryr(l),cp50n(l),cp50s(l),cr50n(l), c * sp50n(l),sp50s(l),sr50n(l) enddo do l=1,3 read(15,'(a5,i2)')c5,ioc c print*,' ioc ',ioc read(15,'(a5)')c5 do n=16,1,-1 read(15,'(f6.1,f5.2,f5.2,f6.0,f6.0,f6.3,f8.1, * f7.1,f6.2,f5.1,f7.2)') * souphi(n),soue(n,ioc),soup(n,ioc),sourr(n,ioc), * sourv(n,ioc),sousp(n,ioc),soua(n,ioc),soudep(n,ioc), * souiep(n,ioc),souir(n,ioc),souiv(n,ioc) c write(6,'(f6.1,f5.2,f5.2,f6.0,f6.0,f6.3,f8.1, c * f7.1,f6.2,f5.1,f7.2)') c * souphi(n),soue(n,ioc),soup(n,ioc),sourr(n,ioc), c * sourv(n,ioc),sousp(n,ioc),soua(n,ioc),soudep(n,ioc), c * souiep(n,ioc),souir(n,ioc),souiv(n,ioc) enddo do n=17,32 read(15,'(f6.1,f5.2,f5.2,f6.0,f6.0,f6.3,f8.1, * f7.1,f6.2,f5.1,f7.2)') * souphi(n),soue(n,ioc),soup(n,ioc),sourr(n,ioc), * sourv(n,ioc),sousp(n,ioc),soua(n,ioc),soudep(n,ioc), * souiep(n,ioc),souir(n,ioc),souiv(n,ioc) c write(6,'(f6.1,f5.2,f5.2,f6.0,f6.0,f6.3,f8.1, c * f7.1,f6.2,f5.1,f7.2)') c * souphi(n),soue(n,ioc),soup(n,ioc),sourr(n,ioc), c * sourv(n,ioc),sousp(n,ioc),soua(n,ioc),soudep(n,ioc), c * souiep(n,ioc),souir(n,ioc),souiv(n,ioc) enddo enddo c c years to which runoff curve should be extrapolated do n=1,neextr yrextr(n)=souryr(newr)+real(n) enddo c c linear/expon. extrapolation of runoff do n=1,neextr c cr50ne(n)=cr50n(newr)+real(n)*(cr50n(newr)-cr50n(newr-1)) cr50ne(n)=cr50n(newr)*exp(-deccon*dt*12.*real(n)) expo=exp(-deccon*dt*12.*real(n)) c print*,' n deccon dt expo ',n,deccon,dt,expo c print*,' n cr50ne cr50n(newr) ',n,cr50ne(n),cr50n(newr) enddo ireawr=1 endif c read Doney's stuff c---read coefficient time series (first two principal c components of doneys's tritium precip. function) c and assoc. spatial patterns c do n=1,nedon read(19,'(f8.1,2(1x,f6.3))')donyr(n), * (cptdon(n,m),m=1,2) c write(6,'(f8.1,2(1x,f6.3))')donyr(n), c * (cptdon(n,m),m=1,2) enddo return end c------------------------------------------------------------ subroutine tracer_init(npt,nz,nzi,nstart,nxp,nyp,iox,trmf,tmf,hmf,xm,ym,trclm) c------------------------------------------------------------ implicit real(a-h,o-z),integer(i-n) include 'comm_para.h' include 'comm_new.h' include 'comm_tracer.h' dimension trmf(npt,nz,1),hmf(npt,nz),tmf(npt,nz),trclm(npt,nz,1) dimension xm(1),ym(1),nzi(npt) dimension iox(1) dimension xga(MAXTR),yga(MAXTR),dga(MAXTR),rga(MAXTR) common /gauss_trac/ xga, yga, dga, rga nptz = npt*nz c If Restart>=1 w/ tracers, c then be sure NOT to reinitialize if (irest.ge.1) return do j=1,ntrac_tot initj = inittr(j) if (initj .eq. 1) then do i=1,npt do k=1,nz trmf(i,k,j) = 0.0 enddo enddo elseif (initj .eq. 2) then do i=1,npt do k=1,nzi(i) trmf(i,k,j) = tr_in(j) enddo enddo elseif (initj .eq. 3 .or. initj.eq.11) then c initialize with data call tr_data_in(j,npt,nz,hmf,trclm(1,1,j)) do i=1,npt do k=1,nzi(i) trmf(i,k,j) = trclm(i,k,j) enddo enddo elseif (initj .eq. 4) then c biological init with temp function call biotr_init(npt,nz,nzi,j-ntrac,trmf(1,1,j),tmf) elseif (initj .eq. 5) then c TR_INIT = latitude call latitude_init(npt,nz,nxp,nyp,iox,trmf(1,1,j),hmf,xm,ym) elseif (initj .ge. 6 .or. initj .le. 8) then call sin_lat_init(npt,nz,nxp,nyp,initj,iox,trmf(1,1,j),hmf,xm,ym) elseif (initj .eq. 10) then c TR_INIT = depth call z_init(npt,nz,trmf(1,1,j),hmf) elseif (initj .eq. 21) then call gauss_init(npt,nz,nxp,nyp, & dga(j),rga(j),yga(j),xga(j),iox,xm,ym,trmf(1,1,j),hmf) endif if (initj .eq. 22) then call lumc_init_a(npt,nz,nxp,nyp,initj,trmf(1,1,j),ym,iox) elseif (initj .eq. 23) then call lumc_init_b(npt,nz,nxp,nyp,initj,trmf(1,1,j),ym,iox) elseif (initj .eq. 24) then call lumc_init_c(npt,nz,nxp,nyp,initj,trmf(1,1,j),ym,iox) endif do i=1,npt do k=nzi(i)+1,nz trmf(i,k,j) = 0.0 enddo enddo enddo return end c-------------------------------------------------- subroutine lumc_init_a(npt,nz,nx,ny,initj,tr,yy,iox) c-------------------------------------------------- implicit real(a-h,o-z),integer(i-n) dimension tr(npt,nz) dimension yy(ny) dimension iox(npt) do i=1,npt do k=1,nz j = (iox(i)-1)/nx + 1 if (yy(j) .le. -18.0) then tr(i,k)=1.0 endif enddo enddo return end c-------------------------------------------------- subroutine lumc_init_b(npt,nz,nx,ny,initj,tr,yy,iox) c-------------------------------------------------- implicit real(a-h,o-z),integer(i-n) dimension tr(npt,nz) dimension yy(ny) dimension iox(npt) do i=1,npt do k=1,nz j = (iox(i)-1)/nx + 1 if ((yy(j) .gt. -18.0).and.(yy(j) .lt. 18.0)) then tr(i,k)=1.0 endif enddo enddo return end c-------------------------------------------------- subroutine lumc_init_c(npt,nz,nx,ny,initj,tr,yy,iox) c-------------------------------------------------- implicit real(a-h,o-z),integer(i-n) dimension tr(npt,nz) dimension yy(ny) dimension iox(npt) do i=1,npt do k=1,nz j = (iox(i)-1)/nx + 1 if (yy(j) .ge. 18.0) then tr(i,k)=1.0 endif enddo enddo return end c ------------------------------------------------------------ subroutine latitude_init(npt,nz,nxp,nyp,iox,trm,hm,xm,ym) c ------------------------------------------------------------ c this subroutine is only appropriate for c a tropical domain (30S <-> 30N) c implicit real(a-h,o-z),integer(i-n) dimension iox(npt) dimension trm(npt,nz),hm(npt,nz) dimension xm(1),ym(1),zz(100) c do i=1,npt j=((iox(i)-1)/nxp)+1 do k=1,nz trm(i,k) = ym(j) enddo enddo return end c ------------------------------------------------------------ subroutine sin_lat_init(npt,nz,nxp,nyp,init_tr,iox,trm,hm,xm,ym) c ------------------------------------------------------------ c this subroutine is only appropriate for c a tropical domain (30S <-> 30N) c implicit real(a-h,o-z),integer(i-n) dimension iox(npt) dimension trm(npt,nz),hm(npt,nz) dimension xm(1),ym(1),zz(100) parameter (RTODEG = 180./3.14159265) c if (init_tr .eq. 6) then do i=1,npt j=((iox(i)-1)/nxp)+1 rlat = ym(j)/rtodeg do k=1,nz if (ym(j).ge. 0.0) then trm(i,k) = sin(rlat) c trm(i,k) = sin(ym(j)) else trm(i,k) = 0.0 endif enddo enddo elseif (init_tr .eq. 7) then do i=1,npt j=((iox(i)-1)/nxp)+1 rlat = ym(j)/rtodeg do k=1,nz if (ym(j).le. 0.0) then trm(i,k) = -sin(rlat) else trm(i,k) = 0.0 endif enddo enddo elseif (init_tr .eq. 8) then do i=1,npt j=((iox(i)-1)/nxp)+1 do k=1,nz trm(i,k) = ym(j) enddo enddo endif return end c ------------------------------------------------------------ subroutine gauss_init(npt,nz,nxp,nyp, & dga,rga,yga,xga,iox,xm,ym,tr,hm) c ------------------------------------------------------------ implicit real(a-h,o-z),integer(i-n) dimension iox(npt), xm(1), ym(1), zz(100), dz(100) dimension tr(npt,nz), hm(npt,nz) do indx=1,npt j = (iox(indx)-1)/nxp + 1 i = iox(indx) - (j-1)*nxp dz(1) = 0.5*hm(indx,1) zz(1) = 0.5*hm(indx,1) termx = exp(-((xm(i)-xga)/rga)**2.) termy = exp(-((ym(j)-yga)/rga)**2.) termz = exp(-zz(1)/dga) tr(indx,1) = termx*termy*termz do k=2,nz dz(k) = hm(indx,k-1) - dz(k-1) zz(k) = zz(k-1) + 2.*dz(k) termz = exp(-zz(k)/dga) tr(indx,k) = termx*termy*termz enddo enddo return end c ------------------------------------------------------------ subroutine z_init(npt,nz,trm,hm) c ------------------------------------------------------------ implicit real(a-h,o-z),integer(i-n) dimension trm(npt,nz),hm(npt,nz) dimension zz(100) c do i=1,npt trm(i,1) = 0.5*hm(i,1) do k=2,nz trm(i,k) = trm(i,k-1) + 0.5*(hm(i,k-1)+hm(i,k)) enddo enddo return end c ------------------------------------------------------- subroutine tr_data_in(itrac,npt,nz,hcm,trclm) c ------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' include 'comm_data.h' include 'comm_tracer.h' real hcm(npt,nz),trclm(npt,nz) nlen = n_tr(itrac) nlen2 = name_tr(itrac) name_temporary = fbtr(itrac) name_temporary2 = ftrnm(itrac) call odb_open(idf_trclim(itrac),name_temporary(1:nlen),0) call odb_rddm(idf_trclim(itrac),'Z',nztr) call odb_rddm(idf_trclim(itrac),'T',ntimes) call mem_alloc(p_ztrclim,nztr,2,'ztrclim') call odb_rdgr(idf_trclim(itrac),'Z',nztr,ztrclim) call data_on_model_grid(idf_trclim(itrac),lclm,name_temporary2(1:nlen2)) call read_linz(idf_trclim(itrac),lclm,npt,mpack,nz,nztr,1,hcm, & trclm, ztrclim,tp,name_temporary2(1:nlen2)) call mem_free(p_ztrclim,nztr,2) return end c ------------------------------------------------------------ subroutine array_init(npt,nz) c ------------------------------------------------------------ include 'comm_new.h' include 'comm_data.h' include 'comm_tracer.h' common /all_loc/ memory_used c c comment: the two arrays tramt(ntrac,nz+1,2) c and trint(ntrac,nz+1,2) c are dimensionalized; c "2" refers to the fact that tracer c amount and tracer variance are calculated; c nptz = npt*nz c nptr = npt*nz*ntrac ncons = 3 nptint = (nz+1)*ntrac*ncons call mem_alloc(p_tramt,nptint,2,'tramt') call mem_alloc(p_trfirst,nptint,2,'trfirst') call mem_alloc(p_vollev,nz,2,'vollev') call mem_alloc(p_trlev,nz,2,'trlev') call mem_alloc(p_sqrlev,nz,2,'sqrlev') return end c ------------------------------------------------------------ subroutine biotr_init(npt,nz,nzi,index,tr,t) c ------------------------------------------------------------ include 'biology.h' CC NPT : number of horizontal grid points CC NZ : number of maximum vertical layers CC NZI : number of vertical layers CC TR : biological tracers array (in mmoleN/m3) CC T : temperature (in Celsius) integer nzi(npt) real tr(npt,nz), t(npt,nz) C C COEFFICIENTS FOR PHYTOPLANKTON FUNCTION C C ZPMAX=0.1 ZPMIN=1.E-2 ZNMIN=6. ZNMAX=1 C ZPMEAN=0.5*(ZPMAX+ZPMIN) ZPDIFF=0.5*(ZPMAX-ZPMIN) ZNMEAN=0.5*(ZNMAX+ZNMIN) ZNDIFF=0.5*(ZNMAX-ZNMIN) ZPISUR2=ASIN(1.) C C RECHERCHE DES MIN ET MAX DE DENSITE : C ZMIN=1.E+10 ZMAX=-1.E+10 DO I=1,NPT DO K=1,NZI(I) ZMIN=MIN(ZMIN,T(I,K)) ZMAX=MAX(ZMAX,T(I,K)) END DO END DO ZMEAN=(ZMIN+ZMAX)*0.5 ZDIFF=(ZMAX-ZMIN)*0.5 if (index .eq. nphy) then DO I=1,NPT DO K=1,NZI(I) TR(I,K)=(ZPMEAN + $ ZPDIFF* SIN(ZPISUR2*(T(I,K)-ZMEAN)/ZDIFF)) END DO END DO elseif (index.eq.nnut) then DO I=1,NPT DO K=1,NZI(I) TR(I,K)=(ZNMEAN + $ ZNDIFF* SIN(ZPISUR2*(T(I,K)-ZMEAN)/ZDIFF)) END DO END DO else print*,'error in biotr_init' stop endif return end