c*********************************************************************** subroutine mzcud( tauinf,tau, c,cup,cdw,vc,taugr, @ ib,isot,icfout,itableout ) c old times mlp first version of mzcf c a.k.murphy method to avoid extrapolation in the curtis matrix c feb-89 malv AKM method to avoid extrapolation in C.M. c 25-sept-96 cristina dejar las matrices en doble precision c jan 98 malv version para mz1d c oct 01 malv update version for fluxes for hb and fb c jul 2011 malv+fgg Adapted to LMD-MGCM c*********************************************************************** implicit none include 'comcstfi.h' include 'nltedefs.h' include 'nlte_atm.h' include 'nlte_data.h' include 'tcr_15um.h' include 'nlte_curtis.h' c arguments real*8 c(nl,nl), cup(nl,nl), cdw(nl,nl) ! o real*8 vc(nl), taugr(nl) ! o real*8 tau(nl,nl) ! i real*8 tauinf(nl) ! i integer ib ! i integer isot ! i integer icfout, itableout ! i c external external bandid character*2 bandid c local variables integer i, in, ir, iw, itblout real*8 cfup(nl,nl), cfdw(nl,nl) real*8 a(nl,nl), cf(nl,nl) character isotcode*2, bcode*2 c formats 101 format(i1) 202 format(i2) 180 format(a80) 181 format(a80) c*********************************************************************** if (isot.eq.1) isotcode = '26' if (isot.eq.2) isotcode = '28' if (isot.eq.3) isotcode = '36' if (isot.eq.4) isotcode = '27' if (isot.eq.5) isotcode = 'co' bcode = bandid( ib ) ! write (*,*) ' ' do in=1,nl do ir=1,nl cf(in,ir) = 0.0d0 cfup(in,ir) = 0.0d0 cfdw(in,ir) = 0.0d0 c(in,ir) = 0.0d0 cup(in,ir) = 0.0d0 cdw(in,ir) = 0.0d0 a(in,ir) = 0.0d0 end do vc(in) = 0.0d0 taugr(in) = 0.0d0 end do c the next lines are a reduced and equivalent way of calculating c the c(in,ir) elements for n=2,nl1 and r=1,nl c do in=2,nl1 c do ir=1,nl c if(ir.eq.1)then c c(in,ir)=tau(in-1,1)-tau(in+1,1) c elseif(ir.eq.nl)then c c(in,ir)=tau(in+1,nl1)-tauinf(in+1)-tau(in-1,nl1)+tauinf(in-1) c else c c(in,ir)=tau(in+1,ir-1)-tau(in+1,ir)-tau(in-1,ir-1)+tau(in-1,ir) c end if c c(in,ir)=c(in,ir)*pi*deltanu(ib)/(2.*deltaz*1.0e5) c end do c end do c go to 1000 c calculation of the matrix cfup(nl,nl) cfup(1,1) = 1.d0 - tau(1,1) do in=2,nl do ir=1,in if (ir.eq.1) then cfup(in,ir) = tau(in,ir) - tau(in,1) elseif (ir.eq.in) then cfup(in,ir) = 1.d0 - tau(in,ir-1) else cfup(in,ir) = tau(in,ir) - tau(in,ir-1) end if end do end do ! contribution to upwards fluxes from bb at bottom : do in=1,nl taugr(in) = tau(in,1) enddo c calculation of the matrix cfdw(nl,nl) cfdw(nl,nl) = 1.d0 - tauinf(nl) do in=1,nl-1 do ir=in,nl if (ir.eq.in) then cfdw(in,ir) = 1.d0 - tau(in,ir) elseif (ir.eq.nl) then cfdw(in,ir) = tau(in,ir-1) - tauinf(in) else cfdw(in,ir) = tau(in,ir-1) - tau(in,ir) end if end do end do c calculation of the matrix cf(nl,nl) do in=1,nl do ir=1,nl if (ir.eq.1) then ! version con l_bb(tg) = l_bb(t(1))=j(1) (see also vc below) ! cf(in,ir) = tau(in,ir) ! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see also vc below) cf(in,ir) = tau(in,ir) - tau(in,1) elseif (ir.eq.nl) then cf(in,ir) = tauinf(in) - tau(in,ir-1) else cf(in,ir) = tau(in,ir) - tau(in,ir-1) end if end do end do c definition of the a(nl,nl) matrix do in=2,nl-1 do ir=1,nl if (ir.eq.in+1) a(in,ir) = -1.d0 if (ir.eq.in-1) a(in,ir) = +1.d0 a(in,ir) = a(in,ir) / ( 2.d0*deltaz*1.d5 ) end do end do ! this is not needed anymore in the akm scheme ! a(1,1) = +3.d0 ! a(1,2) = -4.d0 ! a(1,3) = +1.d0 ! a(nl,nl) = -3.d0 ! a(nl,nl1) = +4.d0 ! a(nl,nl2) = -1.d0 c calculation of the final curtis matrix ("reduced" by murphy's method) if (isot.ne.5) then do in=1,nl do ir=1,nl cf(in,ir) = cf(in,ir) * pi*deltanu(isot,ib) cfup(in,ir) = cfup(in,ir) * pi*deltanu(isot,ib) cfdw(in,ir) = cfdw(in,ir) * pi*deltanu(isot,ib) end do taugr(in) = taugr(in) * pi*deltanu(isot,ib) end do else do in=1,nl do ir=1,nl cf(in,ir) = cf(in,ir) * pi*deltanuco enddo taugr(in) = taugr(in) * pi*deltanuco enddo endif do in=2,nl-1 do ir=1,nl do i=1,nl ! only c contains the matrix a. matrixes cup,cdw dont because ! these two will be used for flux calculations, not ! only for flux divergencies c(in,ir) = c(in,ir) + a(in,i) * cf(i,ir) ! from this matrix we will extract (see below) the ! nl2 x nl2 "core" for the "reduced" final curtis matrix. end do cup(in,ir) = cfup(in,ir) cdw(in,ir) = cfdw(in,ir) end do ! version con l_bb(tg) = l_bb(t(1))=j(1) (see cf above) !vc(in) = c(in,1) ! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see cf above) if (isot.ne.5) then vc(in) = pi*deltanu(isot,ib)/( 2.d0*deltaz*1.d5 ) * @ ( tau(in-1,1) - tau(in+1,1) ) else vc(in) = pi*deltanuco/( 2.d0*deltaz*1.d5 ) * @ ( tau(in-1,1) - tau(in+1,1) ) endif end do 5 continue ! write (*,*) 'mztf/1/ c(2,*) =', (c(2,i), i=1,nl) ! call elimin_dibuja(c,nl,itableout) c ventana del smoothing de c es nw=3 y de vc es 5 (puesto en lisa): c subroutine elimin_mz4(c,vc,ilayer,nl,nan,iw, nw) iw = nan if (isot.eq.4) iw = 5 ! eliminates values < 1.d-19 if (itableout.eq.30) then itblout = 0 else itblout = itableout endif call elimin_mz1d (c,vc,0,iw,itblout,nw) ! upper boundary condition ! j'(nl) = j'(nl1) ==> j(nl) = 2j(nl1) - j(nl2) ==> do in=2,nl-1 c(in,nl-2) = c(in,nl-2) - c(in,nl) c(in,nl-1) = c(in,nl-1) + 2.d0*c(in,nl) cup(in,nl-2) = cup(in,nl-2) - cup(in,nl) cup(in,nl-1) = cup(in,nl-1) + 2.d0*cup(in,nl) cdw(in,nl-2) = cdw(in,nl-2) - cdw(in,nl) cdw(in,nl-1) = cdw(in,nl-1) + 2.d0*cdw(in,nl) end do ! j(nl) = j(nl1) ==> ! do in=2,nl1 ! c(in,nl1) = c(in,nl1) + c(in,nl) ! end do ! 1000 continue if (icfout.eq.1) then ! if (ib.eq.1 .or. ib.eq.12 .or. ib.eq.16 .or. ib.eq.18) then ! codmatrx = codmatrx_fb ! else ! codmatrx = codmatrx_hot ! end if ! if (ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! ibcode2 = '0'//ibcode1 ! else ! write ( ibcode2, 202) ib ! endif ! open ( 1, access='sequential', form='unformatted', file= ! @ dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat') ! open ( 2, access='sequential', form='unformatted', file= ! @ dircurtis//'cflup'//isotcode//dn//ibcode2//codmatrx//'.dat') ! open ( 3, access='sequential', form='unformatted', file= ! @ dircurtis//'cfldw'//isotcode//dn//ibcode2//codmatrx//'.dat') ! open ( 4, access='sequential', form='unformatted', file= ! @ dircurtis//'cflgr'//isotcode//dn//ibcode2//codmatrx//'.dat') ! write(1) dummy ! write(1) ' format: (vc(n),(ch(n,r),r=2,nl-1),n=2,nl-1)' ! do in=2,nl-1 ! write(1) vc(in), (c(in,ir) , ir=2,nl-1 ) !! write (*,*) in, vc(in) ! end do ! write(2) dummy ! write(2)' format: (cfup(n,r),r=1,nl), n=1,nl)' ! do in=1,nl ! write(2) ( cup(in,ir) , ir=1,nl ) ! end do ! write(3) dummy ! write(3)' format: (cfdw(n,r),r=1,nl), n=1,nl)' ! do in=1,nl ! write(3) (cdw(in,ir) , ir=1,nl ) ! end do ! write(4) dummy ! write(4)' format: (taugr(n), n=1,nl)' ! do in=1,nl ! write(4) (taugr(in), ir=1,nl ) ! end do ! !write (*,*) ' Last value in file: ', taugr(nl) ! write (*,'(1x,30hcurtis matrix written out in: ,a,a,a,a)' ) ! @ dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat', ! @ dircurtis//'cflup'//isotcode//dn//ibcode2//codmatrx//'.dat', ! @ dircurtis//'cfldw'//isotcode//dn//ibcode2//codmatrx//'.dat', ! @ dircurtis//'cflgr'//isotcode//dn//ibcode2//codmatrx//'.dat' ! close (1) ! close (2) ! close (3) ! close (4) else ! write (*,*) ' no curtis matrix output file ', char(10) end if if (itableout.eq.30) then ! Force output of C.M. in ascii file ! if (ib.eq.1 .or. ib.eq.12 .or. ib.eq.16 .or. ib.eq.18) then ! codmatrx = codmatrx_fb ! else ! codmatrx = codmatrx_hot ! end if ! if (ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! ibcode2 = '0'//ibcode1 ! else ! write ( ibcode2, 202) ib ! endif ! open (10, file= ! & dircurtis//'table'//isotcode//dn//ibcode2//codmatrx//'.dat') ! write(10,*) nl, ' = number of layers ' ! write(10,*) ' format: (vc(n),(ch(n,r),r=2,nl-1),n=2,nl-1)' ! do in=2,nl-1 ! write(10,*) vc(in), (c(in,ir) , ir=2,nl-1 ) ! enddo ! close (10) endif c end return end