!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! ! ! This program reads ocean and ice data on tripolar grid ! ! and interpolate the data onto regular lat/lon grid. ! ! It contains 5 subroutines: read_ocn_tp, read_ice_tp, tripo_interp40tc, ! ! tripo_interp1tc, tripo ! ! ! ! Xingren Wu (Xingren.Wu@noaa.gov) ! ! May 18, 2007 ! ! Xingren Wu ! ! Nov 16, 2007 - modified ! ! fix kpds(13) to kpds(15) for handle different forecast time unit ! ! Xingren Wu ! ! Dec 4, 2007 - modified ! ! add extra fields to match production output ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! module tripo_intp_mod implicit none contains subroutine read_ocn_tp(im,jm,km,geolon_c,geolat_c,geolon_t,geolat_t, & t,s,u,v,w,eta,sfcflx,pme,mld,taux,tauy,dckt,dcks,vfc,ocean_file) include "netcdf.inc" integer im,jm,km,i,j integer status, ncid integer :: id_geolon_c,id_geolon_t,id_geolat_c,id_geolat_t, & id_temp,id_salt,id_u,id_v,id_wt, & id_eta_t,id_sfc_hflux,id_pme,id_mld,id_tau_x,id_tau_y integer :: id_diff_cbt_kpp_t,id_diff_cbt_kpp_s,id_vert_fric_coeff real, dimension(im,jm,km) :: t,s,u,v,w real, dimension(im,jm,km) :: dckt,dcks,vfc real, dimension(im,jm) :: eta,sfcflx,pme,mld,taux,tauy, & geolon_c,geolon_t,geolat_c,geolat_t character(len=120) :: ocean_file status = nf_open(ocean_file, NF_NOWRITE, ncid) status = nf_inq_varid(ncid, 'temp ', id_temp) status = nf_inq_varid(ncid, 'salt ', id_salt) status = nf_inq_varid(ncid, 'u' , id_u ) status = nf_inq_varid(ncid, 'v' , id_v ) status = nf_inq_varid(ncid, 'wt ' , id_wt) status = nf_inq_varid(ncid, 'eta_t', id_eta_t) status = nf_inq_varid(ncid, 'sfc_hflux', id_sfc_hflux) status = nf_inq_varid(ncid, 'pme ', id_pme) status = nf_inq_varid(ncid, 'mld', id_mld) status = nf_inq_varid(ncid, 'tau_x', id_tau_x) status = nf_inq_varid(ncid, 'tau_y', id_tau_y) status = nf_inq_varid(ncid, 'diff_cbt_kpp_t',id_diff_cbt_kpp_t) status = nf_inq_varid(ncid, 'diff_cbt_kpp_s',id_diff_cbt_kpp_s) status = nf_inq_varid(ncid, 'vert_fric_coeff',id_vert_fric_coeff) status = nf_inq_varid(ncid, 'geolon_c' , id_geolon_c) status = nf_inq_varid(ncid, 'geolon_t' , id_geolon_t) status = nf_inq_varid(ncid, 'geolat_c' , id_geolat_c) status = nf_inq_varid(ncid, 'geolat_t' , id_geolat_t) ! read variable status = nf_get_var_real(ncid, id_geolon_c, geolon_c) status = nf_get_var_real(ncid, id_geolon_t, geolon_t) status = nf_get_var_real(ncid, id_geolat_c, geolat_c) status = nf_get_var_real(ncid, id_geolat_t, geolat_t) status = nf_get_var_real(ncid, id_temp, t) !temp status = nf_get_var_real(ncid, id_salt, s) !salt status = nf_get_var_real(ncid, id_u, u) !u status = nf_get_var_real(ncid, id_v, v) !v status = nf_get_var_real(ncid, id_wt, w) !wt status = nf_get_var_real(ncid, id_eta_t, eta) !eta_t status = nf_get_var_real(ncid, id_sfc_hflux, sfcflx) !sfc_hflux status = nf_get_var_real(ncid, id_pme, pme) !pme status = nf_get_var_real(ncid, id_mld, mld) !mld status = nf_get_var_real(ncid, id_tau_x, taux) !tau_x status = nf_get_var_real(ncid, id_tau_y, tauy) !tau_y status = nf_get_var_real(ncid,id_diff_cbt_kpp_t,dckt) !diff_cbt_kpp_t status = nf_get_var_real(ncid,id_diff_cbt_kpp_s,dcks) !diff_cbt_kpp_s status = nf_get_var_real(ncid,id_vert_fric_coeff,vfc) !vert_fric_coeff return end subroutine read_ocn_tp subroutine read_ice_tp(im,jm,ki,sinrot,cosrot,& hi,hs,ts,t1,t2,fi,alb,ui,vi,sst,saltf,icefile) include "netcdf.inc" character(len=120) :: icefile integer :: status,ncid integer :: im,jm,ki,i,j integer :: id_hi,id_hs,id_ts,id_t1,id_t2,id_cn,id_alb,id_ui,id_vi, & id_sst,id_saltf,id_sinrot,id_cosrot real*8, dimension(im,jm) :: hi,hs,ts,t1,t2,fi,alb,ui,vi,sst,saltf,cosrot,sinrot real*8, dimension(im,jm,ki) :: cn print *, 'icefile:', icefile status = nf_open(icefile, NF_NOWRITE, ncid) print *, 'status:', status if(status.ne.NF_NOERR) stop 'open' status = nf_inq_varid(ncid, 'HI ', id_hi) status = nf_inq_varid(ncid, 'HS ', id_hs) status = nf_inq_varid(ncid, 'TS ', id_ts) status = nf_inq_varid(ncid, 'T1 ', id_t1) status = nf_inq_varid(ncid, 'T2 ', id_t2) status = nf_inq_varid(ncid, 'CN ', id_cn) status = nf_inq_varid(ncid, 'ALB ', id_alb) status = nf_inq_varid(ncid, 'UI ', id_ui) status = nf_inq_varid(ncid, 'VI ', id_vi) status = nf_inq_varid(ncid, 'SST ', id_sst) status = nf_inq_varid(ncid, 'SALTF ', id_saltf) status = nf_inq_varid(ncid, 'SINROT ', id_sinrot) status = nf_inq_varid(ncid, 'COSROT ', id_cosrot) ! ! read hi,hs,ts,t1,t2,alb,ui,vi,sst,saltf ! status = nf_get_var_double(ncid, id_hi ,hi) status = nf_get_var_double(ncid, id_hs ,hs) status = nf_get_var_double(ncid, id_ts ,ts) status = nf_get_var_double(ncid, id_t1 ,t1) status = nf_get_var_double(ncid, id_t2 ,t2) status = nf_get_var_double(ncid, id_cn ,cn) status = nf_get_var_double(ncid, id_alb ,alb) status = nf_get_var_double(ncid, id_ui ,ui) status = nf_get_var_double(ncid, id_vi ,vi) status = nf_get_var_double(ncid, id_sst ,sst) status = nf_get_var_double(ncid, id_saltf ,saltf) status = nf_get_var_double(ncid, id_sinrot ,sinrot) status = nf_get_var_double(ncid, id_cosrot ,cosrot) fi = sum(cn,3) fi(:,:) = fi(:,:)-cn(:,:,1) do j=1,jm do i=1,im if (hi(i,j).lt.1e-5) then ts(i,j)=sst(i,j) t1(i,j)=sst(i,j) t2(i,j)=sst(i,j) endif if (fi(i,j).gt.1.00001) then print *,'Warning: fi>1:',fi(i,j) fi(i,j)=1.0 endif enddo enddo fi = sum(cn,3) return end subroutine read_ice_tp subroutine tripo_interp40tc(im,jm,ko,im1,jm1,im2,jm2,im3,jm3,im4,jm4, & imtp,jmtp,imo,jmo,imos,imoe,imop1,data,dato,tgrid) use mod_interp integer im,jm,ko integer im1,im2,im3,im4,imtp,imo,imos,imoe,imop1 integer jm1,jm2,jm3,jm4,jmtp,jmo real spv parameter(spv=-1.e+34) ! integer i,j,k ! real, dimension(im,jm,ko) :: data real, dimension(imo,jmo,ko) :: dato ! real*8, dimension(im1,jm1) :: lon1,lat1,var1 real*8, dimension(im2,jm2) :: lon2,lat2,var2 real*8, dimension(im3,jm3) :: lon3,lat3,var3 real*8, dimension(im4,jm4) :: lon4,lat4,var4 real*8, dimension(imo,jmo) :: lont,latt,vart real*8, dimension(imop1,jmo) :: lonc,latc,varc real*8, dimension(imtp,jmtp) :: vartp ! integer, dimension(im1,jm1,ko) :: mask1,IG1,JG1 integer, dimension(im2,jm2,ko) :: mask2 integer, dimension(im3,jm3,ko) :: mask3,IG3,JG3 integer, dimension(im4,jm4,ko) :: mask4 integer, dimension(imo,jmo) :: maskt,num_levels integer, dimension(imop1,jmo) :: maskc integer, dimension(imop1,jmo,ko) :: IG,JG integer :: ierr character(len=120) :: ocnintp_grid logical tgrid ! print *,'tgrid=', tgrid if (tgrid) then ocnintp_grid='OCNINTPCOEF.T' else ocnintp_grid='OCNINTPCOEF.C' endif ! open(21,file=ocnintp_grid,status='old',form='unformatted') rewind (21) read (21) lon1 read (21) lat1 read (21) lon2 read (21) lat2 read (21) lon3 read (21) lat3 read (21) lon4 read (21) lat4 read (21) lont read (21) latt read (21) lonc read (21) latc do k=1,ko read (21) mask1(:,:,k) read (21) mask2(:,:,k) read (21) mask3(:,:,k) read (21) mask4(:,:,k) read (21) IG1(:,:,k) read (21) JG1(:,:,k) read (21) IG3(:,:,k) read (21) JG3(:,:,k) read (21) IG(:,:,k) read (21) JG(:,:,k) enddo read (21) num_levels close (21) ! open(50,file='num_levels_1x1.data',access='direct',recl=imo*jmo*4,form='unformatted') ! read(50,rec=1) num_levels ! close (50) do k=1,ko var2=spv var4=spv do j=1,jm1 do i=1,im1 var1(i,j)=data(i,j,k) enddo enddo do j=1,jmtp do i=1,imtp vartp(i,j)=data(i,jm1+j-1,k) enddo enddo do i=1,jm3 do j=1,im3/2 var3(j,i)=vartp(i,j) var3(jmtp+j,i)=vartp(im-i+1,jmtp-j+1) enddo enddo call INTERP(im1,jm1,im2,jm2,lon1,lat1,lon2,lat2,mask1(:,:,k),mask2(:,:,k),IG1(:,:,k),JG1(:,:,k),var1,var2,ierr,NOSPVAL=.true.) call INTERP(im3,jm3,im4,jm4,lon3,lat3,lon4,lat4,mask3(:,:,k),mask4(:,:,k),IG3(:,:,k),JG3(:,:,k),var3,var4,ierr,NOSPVAL=.true.) ! do i=1,imo do j=1,jm2 vart(i,j)=var2(i,j) enddo do j=jm2+1,jmo vart(i,j)=var4(i,j-jm2) enddo enddo if (.NOT. tgrid) then do j=1,jmo do i=2,imop1 maskc(i,j)=0 varc(i,j)=vart(i-1,j) if (varc(i,j).LE.spv) then varc(i,j)=spv maskc(i,j)=1 endif enddo varc(1,j)=varc(imop1,j) maskc(1,j)=maskc(imop1,j) enddo vart=spv maskt=0 call INTERP(imop1,jmo,imo,jmo,lonc,latc,lont,latt,maskc,maskt,IG(:,:,k),JG(:,:,k),varc,vart,ierr,NOSPVAL=.true.) endif ! QC do j=1,jmo do i=1,imo if (num_levels(i,j) .LT. k) vart(i,j)=spv enddo enddo ! ! arrange data so it is for: 0 -> 360, 90 N -> 90S and 4478 -> 0 ! do j=1,jmo do i=1,imos dato(i,j,ko-k+1)=vart(i+imoe,jmo-j+1) enddo do i=imos+1,imo dato(i,j,ko-k+1)=vart(i-imos,jmo-j+1) enddo enddo enddo ! return end subroutine tripo_interp40tc subroutine tripo_interp1tc(im,jm,ko,im1,jm1,im2,jm2,im3,jm3,im4,jm4, & imtp,jmtp,imo,jmo,imos,imoe,imop1,data,dato,tgrid) use mod_interp integer im,jm,ko integer im1,im2,im3,im4,imtp,imo,imos,imoe,imop1 integer jm1,jm2,jm3,jm4,jmtp,jmo real spv parameter(spv=-1.e+34) ! integer i,j,k ! real, dimension(im,jm) :: data real, dimension(imo,jmo) :: dato ! real*8, dimension(im1,jm1) :: lon1,lat1,var1 real*8, dimension(im2,jm2) :: lon2,lat2,var2 real*8, dimension(im3,jm3) :: lon3,lat3,var3 real*8, dimension(im4,jm4) :: lon4,lat4,var4 real*8, dimension(imo,jmo) :: lont,latt,vart real*8, dimension(imop1,jmo) :: lonc,latc,varc real*8, dimension(imtp,jmtp) :: vartp ! integer, dimension(im1,jm1,ko) :: mask1,IG1,JG1 integer, dimension(im2,jm2,ko) :: mask2 integer, dimension(im3,jm3,ko) :: mask3,IG3,JG3 integer, dimension(im4,jm4,ko) :: mask4 integer, dimension(imo,jmo) :: maskt,num_levels integer, dimension(imop1,jmo) :: maskc integer, dimension(imop1,jmo,ko) :: IG,JG integer :: ierr character(len=120) :: ocnintp_grid logical tgrid ! ! print *,'tgrid=', tgrid if (tgrid) then ocnintp_grid='OCNINTPCOEF.T' else ocnintp_grid='OCNINTPCOEF.C' endif ! open(21,file=ocnintp_grid,status='old',form='unformatted') rewind (21) read (21) lon1 read (21) lat1 read (21) lon2 read (21) lat2 read (21) lon3 read (21) lat3 read (21) lon4 read (21) lat4 read (21) lont read (21) latt read (21) lonc read (21) latc do k=1,ko read (21) mask1(:,:,k) read (21) mask2(:,:,k) read (21) mask3(:,:,k) read (21) mask4(:,:,k) read (21) IG1(:,:,k) read (21) JG1(:,:,k) read (21) IG3(:,:,k) read (21) JG3(:,:,k) read (21) IG(:,:,k) read (21) JG(:,:,k) enddo read (21) num_levels close (21) ! open(50,file='num_levels_1x1.data',access='direct',recl=imo*jmo*4,form='unformatted') ! read(50,rec=1) num_levels ! close (50) var2=spv var4=spv do j=1,jm1 do i=1,im1 var1(i,j)=data(i,j) enddo enddo do j=1,jmtp do i=1,imtp vartp(i,j)=data(i,jm1+j-1) enddo enddo do i=1,jm3 do j=1,im3/2 var3(j,i)=vartp(i,j) var3(jmtp+j,i)=vartp(im-i+1,jmtp-j+1) enddo enddo call INTERP(im1,jm1,im2,jm2,lon1,lat1,lon2,lat2,mask1(:,:,1),mask2(:,:,1),IG1(:,:,1),JG1(:,:,1),var1,var2,ierr,NOSPVAL=.true.) call INTERP(im3,jm3,im4,jm4,lon3,lat3,lon4,lat4,mask3(:,:,1),mask4(:,:,1),IG3(:,:,1),JG3(:,:,1),var3,var4,ierr,NOSPVAL=.true.) ! do i=1,imo do j=1,jm2 vart(i,j)=var2(i,j) enddo do j=jm2+1,jmo vart(i,j)=var4(i,j-jm2) enddo enddo if (.NOT. tgrid) then do j=1,jmo do i=2,imop1 maskc(i,j)=0 varc(i,j)=vart(i-1,j) if (varc(i,j).LE.spv) then varc(i,j)=spv maskc(i,j)=1 endif enddo varc(1,j)=varc(imop1,j) maskc(1,j)=maskc(imop1,j) enddo vart=spv maskt=0 call INTERP(imop1,jmo,imo,jmo,lonc,latc,lont,latt,maskc,maskt,IG(:,:,1),JG(:,:,1),varc,vart,ierr,NOSPVAL=.true.) endif ! QC do j=1,jmo do i=1,imo if (num_levels(i,j) .LT. 1) vart(i,j)=spv enddo enddo ! ! arrange data so it is for: 0 -> 360 and 90 N -> 90S ! do j=1,jmo do i=1,imos dato(i,j)=vart(i+imoe,jmo-j+1) enddo do i=imos+1,imo dato(i,j)=vart(i-imos,jmo-j+1) enddo enddo ! return end subroutine tripo_interp1tc subroutine tripo(idbug,im,jm,km,icefile,ocnfile,iyr,imth,iday,ihr, & mfh,mfhout,flonw,flone,dlon,flatn,flats,dlat,jmtp,imos, & imo,jmo,jmtpo,mkmoc,nreg,tripolat,outfile,mocfile,mfcstcpl,igenocnp) integer nvar,ko,ki,kk,ndtc,mkmoc,nreg,mfcstcpl,igenocnp real hfc parameter(nvar=30,ko=40,ki=5) parameter(ndtc=7) ! character*120 icefile,ocnfile,outfile,mocfile ! integer im,jm,km,jmtp integer im1c,im3c,imtpc,jm1c,jm3c,jmtpc integer im1t,im3t,imtpt,jm1t,jm3t,jmtpt integer im2c,im4c,imtpoc,jm2c,jm4c,jmtpoc integer im2t,im4t,imtpot,jm2t,jm4t,jmtpot integer imo,imop1,jmo,jmtpo,imos,imoe integer idbug,iyr,imth,iday,ihr,mfh,mfhout real dlat,dlon,flats,flatn,flonw,flone,tripolat,dtripolat real factor,undef,spv_tau,val integer i,j,k,ierr,ilev,ind,iret,nv,ndata,nr,nundef,kreg,jtripolat ! real*8, dimension(im,jm) :: hi,hs,ts,t1,t2,fi,alb,ui,vi,sst,saltf,crot,srot ! real, dimension(im,jm,km) :: t,s,u,v,w,vv real, dimension(im,jm,km) :: dckt,dcks,vfc real, dimension(im,jm) :: eta,sfcflx,pme,mld,taux,tauy,uice,vice real, dimension(im,jm) :: clat,clon,tlat,tlon,cosrot,sinrot,varice ! real, dimension(imo,jmo) :: grid,varsfc real, dimension(imo,jmo,km) :: varocn,tocn,socn ! real, dimension(jm,km,nreg) :: vmoc real, dimension(im) :: xv real, dimension(jm) :: yv real zt(ko) real zw(ko) real dtc(ndtc) ! real, dimension(nvar) :: fac integer, dimension(nvar) :: kpds5,kpds6,kpds7,kpds22 integer, dimension(ko) :: levs integer, parameter :: kpds_dim=200 integer, dimension(kpds_dim) :: KPDS,KGDS,JPDS,JGDS logical*1 lbms(imo,jmo) logical :: climate = .false. ! data dtc/2.5,5.,10.,15.,20.,25.,28./ ! ! 13=Potential temperature (2) ! 88=Salinity (2) ! 49=u-component of current (2) ! 50=v-component of current (2) ! 40=Geometric Vertical velocity (2) ! 124=Momentum flux, u component (2) ! 125=Momentum flux, v component (2) ! 198=Sea Surface Height Relative to Geoid (129) ! 91=Ice concentration (2) ! 92=Ice thickness (2) ! 66=Snow depth (2) ! 11=Temperature (2) ! 95=u-component of ice drift (2) ! 96=v-component of ice drift (2) ! 188=Evaporation - Precipitation (2) ! 202=Total downward heat flux at surface (downward is positive) (129) ! 195=Geometric Depth Below Sea Surface (129) ! 197=Ocean Heat Content (129) ! 194=Tropical Cyclone Heat Potential (129) ! 195=Ocean Vertical Heat Diffusivity (128) ! 196=Ocean Vertical Salt Diffusivity (128) ! 197=Ocean Vertical Momentum Diffusivity (128) ! data kpds5/ 13, 88, 49, 50, 40,124,125,198, 91, 92, & 66, 11, 95, 96,188,202,195,195,197,194, & 195,195,195,195,195,195,195,195,196,197/ data kpds6/ 160,160,160,160,160, 1, 1, 1, 1, 1, & 1, 1, 1, 1, 1, 1,237,238,236,239, & 235,235,235,235,235,235,235,160,160,160/ data kpds7/ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0, 0, 0, 30,260, & 25, 50,100,150,200,250,280, 0, 0, 0/ data fac/ 273.15,0.001,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0, & 1.0,273.15,1.0,1.0,-8.64e6,1.0,1.0,1.0,1.0,1.0, & 1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0/ data kpds22/ 2, 5, 3, 3, 9, 3, 3, 3, 3, 2, & 3, 2, 3, 3, 3, 3, 0, 0, -5, -4, & 0, 0, 0, 0, 0, 0, 0, 5, 5, 5/ ! data levs/4478,3972,3483,3016,2579,2174,1807,1479,1193,949, & 747, 584, 459, 366, 303, 262, 238, 225, 215,205, & 195, 185, 175, 165, 155, 145, 135, 125, 115,105, & 95, 85, 75, 65, 55, 45, 35, 25, 15, 5/ ! data spv_tau/-1.0E+5/ data undef/-1.0E+34/ ! data zt/ 5., 15., 25., 35., 45., 55., 65., 75., & 85., 95., 105., 115., 125., 135., 145., 155., & 165., 175., 185., 195., 205., 215., 225., 238., & 262., 303., 366., 459., 584., 747., 949.,1193., & 1479.,1807.,2174.,2579.,3016.,3483.,3972.,4478./ ! data zw/ 10., 20., 30., 40., 50., 60., 70., 80., & 90., 100., 110., 120., 130., 140., 150., 160., & 170., 180., 190., 200., 210., 220., 232., 250., & 283., 335., 413., 522., 666., 848.,1072.,1337., & 1643.,1991.,2377.,2798.,3250.,3728.,4225.,4737./ ! im1c=im jm1c=jm-jmtp+1 imtpc=im jmtpc=jmtp im3c=jmtpc*2 jm3c=imtpc/2 ! im1t=im jm1t=jm-jmtp+1 imtpt=im jmtpt=jmtp im3t=jmtpt*2 jm3t=imtpt/2 ! imoe=imo-imos imop1=imo+1 imtpoc=imo jmtpoc=jmtpo imtpot=imo jmtpot=jmtpo ! im2c=imo jm2c=jmo-jmtpoc im4c=im2c jm4c=jmtpoc ! im2t=imo jm2t=jmo-jmtpot im4t=im2t jm4t=jmtpot ! if (mfcstcpl.eq.1) climate = .true. ! print *,'im = ',im print *,'jm = ',jm print *,'km = ',km print *,'jmtp = ',jmtp print *,'im1c = ',im1c print *,'jm1c = ',jm1c print *,'im3c = ',im3c print *,'jm3c = ',jm3c print *,'imtpc = ',imtpc print *,'jmtpc = ',jmtpc print *,'im1t = ',im1t print *,'jm1t = ',jm1t print *,'im3t = ',im3t print *,'jm3t = ',jm3t print *,'imtpt = ',imtpt print *,'jmtpt = ',jmtpt print *,'im2c = ',im2c print *,'jm2c = ',jm2c print *,'im4c = ',im4c print *,'jm4c = ',jm4c print *,'jmtpoc = ',jmtpoc print *,'im2t = ',im2t print *,'jm2t = ',jm2t print *,'im4t = ',im4t print *,'jm4t = ',jm4t print *,'jmtpot = ',jmtpot print *,'imo = ',imo print *,'jmo = ',jmo print *,'imoe = ',imoe print *,'IGEN_OCNP = ',igenocnp print *,'mfcstcpl = ',mfcstcpl print *,'climate = ',climate call read_ocn_tp(im,jm,km,clon,clat,tlon,tlat,t,s,u,v,w, & eta,sfcflx,pme,mld,taux,tauy,dckt,dcks,vfc,ocnfile) print *,'after call read_ocn_tp' call read_ice_tp(im,jm,ki,srot,crot,hi,hs,ts,t1,t2,fi,alb,ui,vi,sst,saltf,icefile) print *,'after call read_ice_tp' cosrot(:,:)=crot(:,:) sinrot(:,:)=srot(:,:) uice(:,:)=ui(:,:) vice(:,:)=vi(:,:) do k=1,km call rotate_data(im,jm,u(:,:,k),v(:,:,k),cosrot,sinrot,undef) enddo call rotate_data(im,jm,taux,tauy,cosrot,sinrot,undef) call rotate_data(im,jm,uice,vice,cosrot,sinrot,undef) if (mkmoc .EQ. 1) then xv(:)=clon(:,1) yv(:)=clat(1,:) jtripolat=jm do j=1,jm if (yv(j) .GT. tripolat) then jtripolat=j go to 191 endif enddo 191 continue dtripolat=yv(jtripolat-1)-yv(jtripolat-2) print *,'jtripolat,dtripolat: ',jtripolat,dtripolat vv(:,:,:)=v(:,:,:) do j=jtripolat,jm yv(j)=yv(j-1)+dtripolat if (yv(j) .GT. 90.) then print *,'WARNING: check yv(j) @', j,yv(j) stop endif vv(:,j,:)=-10. enddo call amoc(vv,xv,yv,zt,im,jm,km,-10.,tripolat,nreg,vmoc) open (71,file=mocfile,form='unformatted') write (71) iyr,imth,iday,ihr,mfh-mfhout,mfh do kreg=1,nreg do k=1,km write (71) vmoc(:,k,kreg) enddo enddo close (71) endif call baopenwt(51,outfile,ierr) if(ierr.ne.0) then print *,'error opening file ',outfile print *,' exit 5' call errexit(5) endif ! kpds(1)=7 if (igenocnp.GT.0) then kpds(2)=igenocnp else kpds(2)=98 endif kpds(3)=10 kpds(4)=192 kpds(8)=mod(iyr-1,100)+1 kpds(9)=imth kpds(10)=iday kpds(11)=ihr kpds(12)=0 if (mfhout == 1) then kpds(13)=1 else if (mfhout == 3) then kpds(13)=10 else if (mfhout == 6) then kpds(13)=11 else if (mfhout == 12) then kpds(13)=12 else if (mfhout == 24) then kpds(13)=2 else print *,'invalid mhout, must be one of (1 3 6 12 24).' stop endif if (climate) then kpds(13)=1 kpds(14)=mfh kpds(15)=0 kpds(16)=10 else if (mfh > 1530) then kpds(13)=1 kpds(14)=mfh kpds(15)=0 kpds(16)=10 else kpds(14)=mfh/mfhout-1 kpds(15)=kpds(14)+1 kpds(16)=3 endif endif kpds(17)=0 kpds(18)=1 kpds(19)=2 kpds(20)=0 kpds(21)=((iyr-1)/100)+1 kpds(23)=4 kpds(24)=0 kpds(25)=32 print*,'kpds:',kpds(1:25) ! kgds(1)=0 kgds(2)=imo kgds(3)=jmo kgds(4)=nint(flatn*1000.) kgds(5)=nint(flonw*1000.) kgds(6)=128 kgds(7)=nint(flats*1000.) kgds(8)=nint(flone*1000.) kgds(9)=nint(dlon*1000.) kgds(10)=nint(dlat*1000.) kgds(11)=0 kgds(12)=0 kgds(13)=0 kgds(14)=0 kgds(15)=0 kgds(16)=0 kgds(17)=0 kgds(18)=0 kgds(19)=0 kgds(20)=255 kgds(21)=0 kgds(22)=0 ! ndata=imo*jmo ! ind=0 do nv=1,5 factor=fac(nv) kpds(22)=kpds22(nv) print *,nv,' factor ',factor,' kpds22 ',kpds(22) !t if (nv.eq.1) then call tripo_interp40tc(im,jm,km,im1t,jm1t,im2t,jm2t,im3t,jm3t,im4t,jm4t, & imtpt,jmtpt,imo,jmo,imos,imoe,imop1,t,varocn,.true.) do k=1,km do j=1,jmo do i=1,imo if (varocn(i,j,k).LE.undef) then tocn(i,j,km-k+1)=undef else tocn(i,j,km-k+1)=varocn(i,j,k)+273.15 endif enddo enddo enddo ! print *,'tocn(360:365,180:185,1)',tocn(360:365,180:185,1) endif !s if (nv.eq.2) then call tripo_interp40tc(im,jm,km,im1t,jm1t,im2t,jm2t,im3t,jm3t,im4t,jm4t, & imtpt,jmtpt,imo,jmo,imos,imoe,imop1,s,varocn,.true.) do k=1,km do j=1,jmo do i=1,imo socn(i,j,k)=varocn(i,j,km-k+1) enddo enddo enddo ! print *,'socn(360:365,180:185,1)',socn(360:365,180:185,1) endif !u if (nv.eq.3) & call tripo_interp40tc(im,jm,km,im1c,jm1c,im2c,jm2c,im3c,jm3c,im4c,jm4c, & imtpc,jmtpc,imo,jmo,imos,imoe,imop1,u,varocn,.false.) !v if (nv.eq.4) & call tripo_interp40tc(im,jm,km,im1c,jm1c,im2c,jm2c,im3c,jm3c,im4c,jm4c, & imtpc,jmtpc,imo,jmo,imos,imoe,imop1,v,varocn,.false.) !w if (nv.eq.5) then call tripo_interp40tc(im,jm,km,im1t,jm1t,im2t,jm2t,im3t,jm3t,im4t,jm4t, & imtpt,jmtpt,imo,jmo,imos,imoe,imop1,w,varocn,.true.) do k=1,km-1 kk=km-k+1 do j=1,jmo do i=1,imo if (varocn(i,j,k).LE.undef) then varocn(i,j,k)=varocn(i,j,k+1) else varocn(i,j,k)=(varocn(i,j,k+1)*(zw(kk)-zt(kk)) & +varocn(i,j,k)*(zt(kk)-zw(kk-1)))/(zw(kk)-zw(kk-1)) endif enddo enddo enddo endif do k=1,km ind=ind+1 ilev=levs(k) ! make bit-map.... do j=1,jmo do i=1,imo val=varocn(i,j,k) if (val.eq.undef) then lbms(i,j) = .false. else if (nv.eq.1) then grid(i,j)=val+factor else grid(i,j)=val*factor endif lbms(i,j) = .true. endif enddo enddo print *,' record written ',ind,nv,k,grid(92,125),lbms(92,125) kpds(5)=kpds5(nv) kpds(6)=kpds6(nv) kpds(7)=ilev CALL PUTGB(51,NDATA,KPDS,KGDS,LBMS,GRID,IRET) if (iret.ne.0) then print *,' error in PUTGB for iret ',iret,(kpds(nr),nr=5,7) print *,' exit 6' call errexit(6) endif ! !.. end level-loop enddo !.. end variable-loop enddo !... now read and grib 5 surface records... ! do nv=6,27 factor=fac(nv) kpds(22)=kpds22(nv) print *,nv,' factor ',factor,' kpds22 ',kpds(22) ind=ind+1 kpds(5)=kpds5(nv) kpds(6)=kpds6(nv) kpds(7)=kpds7(nv) if (nv .EQ. 8 .or. (nv.GE.16 .AND. nv.LE.27)) then kpds(19)=129 else if (nv .EQ. 15) then kpds(19)=128 else kpds(19)=2 endif ! taux if (nv.eq.6) then call tripo_interp1tc(im,jm,km,im1c,jm1c,im2c,jm2c,im3c,jm3c,im4c,jm4c, & imtpc,jmtpc,imo,jmo,imos,imoe,imop1,taux,varsfc,.false.) do j=1,jmo do i=1,imo if (varsfc(i,j) .LE. spv_tau) varsfc(i,j)=undef enddo enddo endif ! tauy if (nv.eq.7) then call tripo_interp1tc(im,jm,km,im1c,jm1c,im2c,jm2c,im3c,jm3c,im4c,jm4c, & imtpc,jmtpc,imo,jmo,imos,imoe,imop1,tauy,varsfc,.false.) do j=1,jmo do i=1,imo if (varsfc(i,j) .LE. spv_tau) varsfc(i,j)=undef enddo enddo endif ! eta if (nv.eq.8) then call tripo_interp1tc(im,jm,km,im1t,jm1t,im2t,jm2t,im3t,jm3t,im4t,jm4t, & imtpt,jmtpt,imo,jmo,imos,imoe,imop1,eta,varsfc,.true.) endif ! fi if (nv.eq.9) then varice=fi call tripo_interp1tc(im,jm,km,im1t,jm1t,im2t,jm2t,im3t,jm3t,im4t,jm4t, & imtpt,jmtpt,imo,jmo,imos,imoe,imop1,varice,varsfc,.true.) endif ! hi if (nv.eq.10) then varice=hi call tripo_interp1tc(im,jm,km,im1t,jm1t,im2t,jm2t,im3t,jm3t,im4t,jm4t, & imtpt,jmtpt,imo,jmo,imos,imoe,imop1,varice,varsfc,.true.) endif ! hs if (nv.eq.11) then varice=hs call tripo_interp1tc(im,jm,km,im1t,jm1t,im2t,jm2t,im3t,jm3t,im4t,jm4t, & imtpt,jmtpt,imo,jmo,imos,imoe,imop1,varice,varsfc,.true.) endif ! ts if (nv.eq.12) then varice=ts call tripo_interp1tc(im,jm,km,im1t,jm1t,im2t,jm2t,im3t,jm3t,im4t,jm4t, & imtpt,jmtpt,imo,jmo,imos,imoe,imop1,varice,varsfc,.true.) endif ! uice if (nv.eq.13) then varice=uice call tripo_interp1tc(im,jm,km,im1c,jm1c,im2c,jm2c,im3c,jm3c,im4c,jm4c, & imtpc,jmtpc,imo,jmo,imos,imoe,imop1,varice,varsfc,.false.) endif ! vice if (nv.eq.14) then varice=vice call tripo_interp1tc(im,jm,km,im1c,jm1c,im2c,jm2c,im3c,jm3c,im4c,jm4c, & imtpc,jmtpc,imo,jmo,imos,imoe,imop1,varice,varsfc,.false.) endif ! pme if (nv.eq.15) then call tripo_interp1tc(im,jm,km,im1t,jm1t,im2t,jm2t,im3t,jm3t,im4t,jm4t, & imtpt,jmtpt,imo,jmo,imos,imoe,imop1,pme,varsfc,.true.) endif ! sfc_flx if (nv.eq.16) then call tripo_interp1tc(im,jm,km,im1t,jm1t,im2t,jm2t,im3t,jm3t,im4t,jm4t, & imtpt,jmtpt,imo,jmo,imos,imoe,imop1,sfcflx,varsfc,.true.) endif ! mld if (nv.eq.17) then call mixed_layer(imo,jmo,km,tocn,socn,zt,varsfc,undef) endif ! sfc isothm layer depth if (nv.eq.18) then call sfc_isothm_layer(imo,jmo,km,tocn,zt,varsfc,undef) endif ! ocean heat content if (nv.eq.19) then call ocean_heat(imo,jmo,km,tocn,socn,zw,zt,varsfc,undef) endif ! tropical cyc heat potential if (nv.eq.20) then call tchp26(imo,jmo,km,tocn,socn,zw,zt,varsfc,undef) endif ! depth of 7 different isotherms... if (nv.ge.21 .AND. nv.le.27) then i=nv-20 call isothm_layer(imo,jmo,km,dtc(i),tocn,zt,varsfc,undef) endif nundef=0 do j=1,jmo do i=1,imo val=varsfc(i,j) if (val.eq.undef) then lbms(i,j) = .false. nundef=nundef+1 else if (nv.EQ.12) then grid(i,j)=val+factor else grid(i,j)=val*factor endif lbms(i,j) = .true. endif enddo enddo if(nundef.eq.imo*jmo) then print *,' record not written because of all undef ', & ind,kpds(5),kpds(6),kpds(7) else print *,' record written ',ind,kpds(5),kpds(6),kpds(7), & grid(92,125),lbms(92,125),factor CALL PUTGB(51,NDATA,KPDS,KGDS,LBMS,GRID,IRET) if(iret.ne.0) then print *,' error in PUTGB for iret ',iret,(kpds(k),k=5,7) print *,' exit 7' call errexit(7) endif endif enddo if (mkmoc .EQ. 1) then do nv=28,nvar factor=fac(nv) kpds(22)=kpds22(nv) print *,nv,' factor ',factor,' kpds22 ',kpds(22) !diff_cbt_kpp_t if (nv.eq.28) & call tripo_interp40tc(im,jm,km,im1t,jm1t,im2t,jm2t,im3t,jm3t,im4t,jm4t, & imtpt,jmtpt,imo,jmo,imos,imoe,imop1,dckt,varocn,.true.) !diff_cbt_kpp_s if (nv.eq.29) & call tripo_interp40tc(im,jm,km,im1t,jm1t,im2t,jm2t,im3t,jm3t,im4t,jm4t, & imtpt,jmtpt,imo,jmo,imos,imoe,imop1,dcks,varocn,.true.) !vert_fric_coeff if (nv.eq.30) & call tripo_interp40tc(im,jm,km,im1c,jm1c,im2c,jm2c,im3c,jm3c,im4c,jm4c, & imtpc,jmtpc,imo,jmo,imos,imoe,imop1,vfc,varocn,.false.) do k=1,km ind=ind+1 ilev=levs(k) ! make bit-map.... do j=1,jmo do i=1,imo val=varocn(i,j,k) if (val.eq.undef) then lbms(i,j) = .false. else if (nv.eq.1) then grid(i,j)=val+factor else grid(i,j)=val*factor endif lbms(i,j) = .true. endif enddo enddo print *,' record written ',ind,nv,k,grid(92,125),lbms(92,125) kpds(5)=kpds5(nv) kpds(6)=kpds6(nv) kpds(7)=ilev kpds(19)=128 CALL PUTGB(51,NDATA,KPDS,KGDS,LBMS,GRID,IRET) if (iret.ne.0) then print *,' error in PUTGB for iret ',iret,(kpds(nr),nr=5,7) print *,' exit 6' call errexit(6) endif ! !.. end level-loop enddo !.. end variable-loop enddo endif ! (mkmoc .EQ. 1) return end subroutine tripo subroutine rotate_data(im,jm,udata,vdata,rotcos,rotsin,undef) integer im,jm,i,j real, dimension(im,jm) :: udata,vdata,rotcos,rotsin real :: u,v,undef do j=1,jm do i=1,im if (udata(i,j) /= undef) then u=udata(i,j)*rotcos(i,j)-vdata(i,j)*rotsin(i,j) v=udata(i,j)*rotsin(i,j)+vdata(i,j)*rotcos(i,j) udata(i,j)=u vdata(i,j)=v endif enddo enddo return end subroutine rotate_data subroutine isothm_layer(im,jm,km,dtc,temp,zlev,zisothm,undef) real, parameter :: c2k=273.15 integer inumc,im,jm,km integer i,j,k real dtc real, dimension(km) :: tz,zlev real, dimension(im,jm) :: zisothm real, dimension(im,jm,km) :: temp real a,b,tc,undef tc=c2k+dtc do j=1,jm do i=1,im zisothm(i,j)=undef if (temp(i,j,1) .GE. tc) then do k=1,km tz(k)=temp(i,j,k) enddo do k=2,km if (tz(k) .LT. -3.0) go to 111 if (tz(k) .LT. tc) then a = (tz(k)-tc) / (tz(k)-tz(k-1)) b = (tc-tz(k-1)) / (tz(k)-tz(k-1)) zisothm(i,j)=a*zlev(k-1)+b*zlev(k) go to 111 endif enddo endif 111 continue enddo enddo return end subroutine isothm_layer subroutine mixed_layer(im,jm,km,temp,salt,zlev,mld,undef) real, parameter :: disot=0.8, c2k=273.15 integer im,jm,km integer i,j,k,kmsk,kbm,kbp,krf real, dimension(km) :: zlev,plev,sa,ta,th,rho real, dimension(im,jm) :: mld real, dimension(im,jm,km) :: temp,salt real a,b,deltarho,dr,rb,undef do k=1,km plev(k) = press(zlev(k),980.0) enddo do j=1,jm do i=1,im kmsk = 0 do k=1,km if (temp(i,j,k).GT.0.0 .AND. salt(i,j,k).GT.0.0) then ta(k) = temp(i,j,k)-c2k sa(k) = salt(i,j,k) kmsk = k endif enddo if (kmsk.EQ.0 .OR. ta(1).LT.-3.0) then mld(i,j)=undef else deltarho = (density(0.0,ta(1)-disot,sa(1)) & - density(0.0,ta(1),sa(1))) do k=1,kmsk th(k) = theta(plev(k),ta(k),sa(k),0.0) rho(k) = density(0.0,th(k),sa(k)) - 1000.0 enddo krf = 1 kbm = 0 kbp = 0 do k = krf,kmsk if ((rho(k)-rho(krf)) .GE. deltarho) then kbp = k exit endif enddo if (kbp .LE. 1) then mld(i,j) = undef else kbm = kbp - 1 rb = rho(krf) + deltarho dr = rho(kbp) - rho(kbm) a = (rho(kbp) - rb) / dr b = (rb - rho(kbm)) / dr mld(i,j) = zlev(kbm)*a + zlev(kbp)*b endif endif enddo enddo end subroutine mixed_layer subroutine sfc_isothm_layer(im,jm,km,temp,zlev,sitd,undef) real, parameter :: disot=0.8 integer im,jm,km integer i,j,k real, dimension(im,jm) :: sitd real, dimension(im,jm,km) :: temp real, dimension(km) :: zlev,tz real a,b,tc,undef do j=1,jm do i=1,im sitd(i,j)=undef if (temp(i,j,1).GE.0.0) then tc=temp(i,j,1)-disot do k=1,km tz(k)=temp(i,j,k) enddo do k=2,km if (tz(k).LT.0.0) go to 112 if (tz(k).LT.tc) then a = (tz(k)-tc) / (tz(k)-tz(k-1)) b = (tc-tz(k-1)) / (tz(k)-tz(k-1)) sitd(i,j) = a*zlev(k-1) + b*zlev(k) go to 112 endif enddo endif 112 continue enddo enddo return end subroutine sfc_isothm_layer subroutine ocean_heat(im,jm,km,temp,salt,zblev,zlev,ocnhc,undef) integer, parameter :: kmh=26 real, parameter :: c2k=273.15 integer im,jm,km integer i,j,k real, dimension(km) :: zblev,zlev,plev real, dimension(im,jm) :: ocnhc real, dimension(im,jm,km) :: salt,temp real dptlyr,rk,sk,tk,undef real pk,rhm,rhp,tempk,zk k=kmh zk=0.5*(300.0+zblev(k-1)) pk=press(zk,980.0) rhm = (zk-zlev(k-1))/(zlev(k)-zlev(k-1)) rhp = (zlev(k)-zk)/(zlev(k)-zlev(k-1)) do k=1,km plev(k) = press(zlev(k),980.0) enddo do j=1,jm do i=1,im ocnhc(i,j)=undef if (temp(i,j,kmh).GT.0.0 .AND. salt(i,j,kmh).GT.0.0) then ocnhc(i,j)=0. do k=1,kmh-1 tk=temp(i,j,k)-c2k sk=salt(i,j,k) rk=density(plev(k),tk,sk) if (k .eq. 1) then dptlyr=zblev(k) else dptlyr=zblev(k)-zblev(k-1) endif ocnhc(i,j)=ocnhc(i,j) + rk*temp(i,j,k)*dptlyr enddo k=kmh tempk=rhp*temp(i,j,k-1) + rhm*temp(i,j,k) tk=tempk - c2k sk=rhp*salt(i,j,k-1) + rhm*salt(i,j,k) rk=density(pk,tk,sk) dptlyr=300.0-zblev(k-1) ocnhc(i,j)=ocnhc(i,j)+rk*tempk*dptlyr ocnhc(i,j)=ocnhc(i,j)*3996. endif enddo enddo return end subroutine ocean_heat subroutine tchp26(im,jm,km,temp,salt,zblev,zlev,ocnhcp,undef) real, parameter :: c2k=273.15, t26=26.0 integer im,jm,km integer i,j,k,k26 real, dimension(km) :: zblev,zlev,plev,tz real, dimension(im,jm) :: ocnhcp,z26isothm real, dimension(im,jm,km) :: salt,temp real dptlyr,rk,sk,tk,undef real rhm,rhp,pk,zk real a,b,skk,skm,tc,z26 logical*1 lbms(im,jm) tc=c2k+t26 do k=1,km plev(k) = press(zlev(k),980.0) enddo do j=1,jm do i=1,im z26isothm(i,j)=undef lbms(i,j)=.false. if (temp(i,j,1) .GE. tc) then do k=1,km tz(k)=temp(i,j,k) enddo k = 1 do while (tz(k).GE.tc) k26 = k if (k.EQ.km) exit k = k + 1 enddo k = k26 if (tz(k) .GT. tc) then if (k.LT.km .AND. tz(k+1).GT.0.0) then k = k + 1 a = (tz(k)-tc) / (tz(k)-tz(k-1)) b = (tc-tz(k-1)) / (tz(k)-tz(k-1)) z26isothm(i,j) = a*zlev(k-1) + b*zlev(k) lbms(i,j)=.true. else if (k.GE.2 .AND. tz(k).LT.tz(k-1)) then a = (tz(k)-tc) / (tz(k)-tz(k-1)) b = (tc-tz(k-1)) / (tz(k)-tz(k-1)) z26 = a*zlev(k-1) + b*zlev(k) if (z26.LE.zblev(k)) then z26isothm(i,j) = z26 lbms(i,j)=.true. endif endif else if (tz(k).EQ.tc) then z26isothm(i,j) = zlev(k) lbms(i,j)=.true. endif endif enddo enddo ! !---------- get ocean heat potential relative to 26C (TCHP) ------------ ! do j=1,jm do i=1,im if (temp(i,j,1) .GT. 0.0) then ocnhcp(i,j)=0.0 else ocnhcp(i,j)=undef cycle endif if (lbms(i,j)) then ! we have water above 26c z26 = z26isothm(i,j) ! ! case where Z26 is within the topmost layer ! if (z26 .LE. zblev(1)) then tk=temp(i,j,1)-c2k if (salt(i,j,1) .GT. 0.0) then sk=salt(i,j,1) else sk=35. ! fake salinity endif rk=density(plev(1),tk,sk) dptlyr=z26 ocnhcp(i,j) = rk*(tk-t26)*dptlyr*3996. ! ! case where z26 is below the top layer and above the bottom ! else k26 = 1 do k=2,km if (z26.GT.zblev(k-1) .AND. z26.LE.zblev(k)) k26=k enddo ocnhcp(i,j)=0.0 do k=1,k26-1 tk=temp(i,j,k)-c2k if (salt(i,j,K) .GT. 0.0) then sk=salt(i,j,k) else sk=35. ! fake salinity endif rk=density(plev(k),tk,sk) if (k .EQ. 1) then dptlyr=zblev(1) else dptlyr=zblev(k)-zblev(k-1) endif ocnhcp(i,j)=ocnhcp(i,j)+rk*(tk-26.0)*dptlyr enddo k=k26 zk=0.5*(z26+zblev(k-1)) pk=press(zk,980.0) rhm = (zk-zlev(k-1))/(zlev(k)-zlev(k-1)) rhp = (zlev(k)-zk)/(zlev(k)-zlev(k-1)) tk=rhp*temp(i,j,k-1) + rhm*temp(i,j,k) - c2k if (salt(i,j,k-1) .GT. 0.0) then skm=salt(i,j,k-1) else skm=35. ! fake salinity endif if (salt(i,j,k) .GT. 0.0) then skk=salt(i,j,k) else skk=35. ! fake salinity endif sk=(rhp*skm + rhm*skk) rk=density(pk,tk,sk) dptlyr=z26-zblev(k-1) ocnhcp(i,j)=ocnhcp(i,j)+rk*(tk-26.0)*dptlyr ocnhcp(i,j)=ocnhcp(i,j)*3996. endif ! ! case where temperature is above 26C down to the bottom ! else if ((temp(i,j,1)-c2k) .GT. t26) then ocnhcp(i,j)=0.0 do k=1,km if (temp(i,j,k) .GT. undef) then tk=temp(i,j,k)-c2k if (salt(i,j,k) .GT. 0.0) then sk=salt(i,j,k) else sk=35. ! fake salinity endif rk=density(plev(k),tk,sk) if (k .EQ. 1) then dptlyr=zblev(1) else dptlyr=zblev(k)-zblev(k-1) endif ocnhcp(i,j)=ocnhcp(i,j)+rk*(tk-26.0)*dptlyr endif enddo ocnhcp(i,j)=ocnhcp(i,j)*3996. endif enddo enddo return end subroutine tchp26 function press(z, g) ! copy from cfs_ocean_time.f and modified ! depth (z) in meters and grav acc'l (g) in cm/sec**2 integer, parameter :: itr=20 integer i real p, a0, z, g, press real(kind=8) :: e, ae, es ! p = z*(1.0076+z*(2.3487e-6-z*1.2887e-11)) e = zeta(p,g)-z ae = abs(e) es = ae*2. do i = 1,itr a0 = 0.972643+p*(1.32696e-5-p*(6.228e-12+p*1.885e-16)) a0 = a0/(1.0+1.83e-5*p) p = p-((g+1.113e-4*p)/a0)*e*0.001 es = ae e = zeta(p,g)-z ae = abs(e) if (ae .le. 0.01) exit enddo ! press = p ! end function press function zeta(p, glat) ! ! copy from cfs_ocean_time.f and modified real p, glat, z, zeta z = ((-3.434e-12*p+1.113e-7)*p+0.712953)*p+14190.7*log(1.0+1.83e-5*p) z = (z/(glat+1.113e-4*p))*1000. zeta = z ! end function zeta function density(prs, tmp, sal) ! ! copy from cfs_ocean_time.f and modified ! Density is in units of kg/m**3 (1 g/cm**3 = 1000 kg/m**3) real density, prs, tmp, sal real p, t, s, kstp, k0, kw, d0, dw ! s = sal t = tmp p = prs/10.00 ! kw = 19652.21+(148.4206-(2.327105-(1.360477e-2-5.155288e-5*t)*t)*t)*t ! k0 = kw+s*(54.6746-(0.603459-(1.09987e-2-6.1670e-5*t)*t)*t) & +sqrt(s*s*s)*(7.944e-2+(1.6483e-2-5.3009e-4*t)*t) ! kstp = k0+p*((3.239908+(1.43713e-3+(1.16092e-4-5.77905e-7*t)*t)*t) & +s*(2.2838e-3-(1.0981e-5+1.6078e-6*t)*t) & +sqrt(s*s*s)*1.91075e-4 & +p*((8.50935e-5-(6.12293e-6-5.2787e-8*t)*t) & -s*(9.9348e-7-(2.0816e-8+9.1697e-10*t)*t))) ! dw = 999.842594+(6.793952e-2-(9.095290e-3-(1.001685e-4 & -(1.120083e-6-6.536332e-9*t)*t)*t)*t)*t ! d0 = dw+s*(0.824493-(4.0899e-3-(7.6438e-5-(8.2467e-7 & -5.3875e-9*t)*t)*t)*t) & -sqrt(s*s*s)*(5.72466e-3-(1.0227e-4-1.6546e-6*t)*t) & +s*s*4.8314e-4 ! density = d0/(1.0-p/kstp) end function density function theta(p, t, s, pref) real(kind=8), parameter :: sqrt2 = 0.7071067811865475 real theta, p,t, s, pref real del_p, del_t1, del_t2, del_t3, del_t4, tp, th del_p = pref-p del_t1 = del_p*atg(p,t,s) tp = t+0.5*del_t1 del_t2 = del_p*atg((p+0.5*del_p),tp,s) tp = t+(-0.5+sqrt2)*del_t1+(1.0-sqrt2)*del_t2 del_t3 = del_p*atg((p+0.5*del_p),tp,s) tp = t-sqrt2*del_t2+(1.0+sqrt2)*del_t3 del_t4 = del_p*atg(pref,tp,s) th = (t+(del_t1+(1.0-sqrt2)*del_t2*2.0 & + (1.0+sqrt2)*del_t3*2.0+del_t4)/6.0) theta = th end function theta function atg(p, t, s) real atg, p, t, s, ds, a ds = s-35.0 a = (((-2.1687e-16*t+1.8676e-14)*t-4.6206e-13)*p & +((2.7759e-12*t-1.1351e-10)*ds+((-5.4481e-14*t & +8.733e-12)*t-6.7795e-10)*t+1.8741e-8))*p & +(-4.2393e-8*t+1.8932e-6)*ds & +((6.6228e-10*t-6.836e-8)*t+8.5258e-6)*t+3.5803e-5 atg = a end function atg end module tripo_intp_mod