module mod_common_ice use mod_xc implicit none ! --- common blocks for the ice and snow part of the model common /icesnw1/ & & albi_m &! -- max albedo over melting ice &,albi_d &! -- max albedo over dry ice &,albi &! -- albedo over ice &,albs &! -- albedo over snow &,albs_m &! -- albedo over melting snow &,albw_d &! -- albedo for diff light over water &,albw &! -- albedo for direct light over water &,rhoice &! kg / m^3 density of ice &,rhosnw &! kg / m^3 density of snow &,rhowat &! kg / m^3 density of pure water &,emiss &! -- emissivity of water &,gasconst &! Pa m^3 / (K kg) gas constant &,rkice &! W / (m K) ice conductivity &,hocond &! J / kg heat of condensation/vap &,fusi &! J / m^3 heat of fusion of ice &,fuss &! J / m^3 heat of fusion of snow !TW changed units -see icedat routine below !-if J/kg, fuss~fusi since both water &,fice_max &! -- maximum fractional ice cover &,tice_m &! K melting point of ice &,tsnw_m &! K melting point of snow &,hice_min &! m minimum ice thickness &,epsmol &! -- molecular weight of h2o/dry air &,albsa ! -- array containing snow albedo save/icesnw1/ common /icesnw2/ & & sice &! per mil salinity of seaice &,rksnw &! W / (m K) snow conductivity &,cpair &! J / (K kg) specific heat of dry air &,cpsw &! J / (K kg) specific heat of seawater &,stefanb &! W / (m^2 K^4) stefan-boltzman constant &,aice &! -- vapor pressure parameters &,bice &! K .. &,awater &! -- .. &,bwater &! K .. &,t0deg ! K zero deg celcius in K save/icesnw2/ real albi_m,albi,albs,albs_m,albw,albw_d, albi_d, & & rhoice,rhosnw,rhowat,emiss,gasconst,rkice,hocond, & & fusi,fuss,fice_max,tice_m,tsnw_m,hice_min,epsmol,sice, & & rksnw,cpair,cpsw,stefanb,aice,bice,awater,bwater,t0deg & & ,albsa(12) real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & & ficem, & & hicem, & & tsrfm, & & hsnwm, & & ticem, & & qbrine, & & cawdir, & & pemnp, & & clat, & & radfl0, & & tauxice, & & tauyice, & & qfrz, & & iceu, & & icev, & & delta_icevol, & & delta_snwvol, & & surf_albedo_sum, & & surf_qsw_sum #if defined (THERM_DIAG) & ,delta_ficem & & ,delta_hicem #endif !TW: add delta_fice,delta_hice to give diagnostics about lat/vert growth/melt real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & MIZ_MASK ! NB if WAVES, dfloe,dfloe_miz, ! dfloe_pack_thresh,dfloe_land ! are now in mod_common_wavesice.F #if defined(ICE_DYN_DIAG) real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & & stressp, & & stressm, & & stress12, & & pice, & & strainI, & & strainII common/iceparam/ stressp,stressm,stress12, & & pice,strainI,strainII save/iceparam/ #if defined (ICE_DYN_DIAG_DUMP) !! option to periodically dump ice diagnostic variables !! to netcdf files real*8, parameter :: icedump_freq = 3.0 !!time between dumps (in hours) !logical, parameter :: dump_time0 = .false. !!dump ice stuff on 1st day? (check restart) logical, parameter :: dump_time0 = .true. !!dump ice stuff on 1st day? (check restart) logical :: DO_INIT_ICEDUMP !!1st call to ice_dyn_diag_dump or not #endif !! sometimes common_blocks.h can't be used to get itest,jtest !! -then you can get these variables by: !! use mod_common_ice, only: mnproc_2,itest_2,jtest_2 !! if (mnproc==mnproc_2) print*,vbl(itest_2,jtest_2) !! NB set them in iceinit below integer :: mnproc_2,itest_2,jtest_2 #endif common/iceparam/ ficem,hicem,tsrfm,hsnwm,ticem, & & qbrine,cawdir,pemnp,clat,radfl0,tauxice,tauyice, & & qfrz,iceU,iceV,delta_icevol,delta_snwvol, & & surf_albedo_sum,surf_qsw_sum save/iceparam/ #if defined (ICE_NEST) && defined (ICE) real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,2) :: & & uicenest, & & vicenest, & & hicenest, & & ficenest, & hsnwnest #endif #if defined (TEST_ICE_AGE) real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & & fy_frac , & & fy_age, & & rdg_frac #endif #if defined (PARAM_EST) real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & & msshb , & & sstb #endif #if defined (ALBSNW_EVOL) real albs_min, albs_max, hsnw_lim real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & & albsnwm, & & surf_fscov real taua, & & tauf, & & swe_newalb #endif #if defined (SSNOWD_ICE) real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & & hmelt ! Cumulated melt depth (m) real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & & hprcp ! Cumulated precipitations depth (m) real cv_snw ! coefficient of variation : depends on the surface. ! Based on SHEBA dataset for sea ice. #endif contains subroutine icedat implicit none albi_m =.60 ! -- max albedo over melting ice albi_d =.73 ! -- max albedo over dry ice albs =.80 ! -- albedo over snow albs_m =.71 ! -- albedo over melting snow albw_d =.065 ! -- albedo over water; diff light rhoice =900. ! kg / m^3 density of ice rhosnw =330. ! kg / m^3 density of snow rhowat =1000. ! kg / m^3 density of pure water emiss =.97 ! -- emissivity of water gasconst =.287e3 ! Pa m^3 / (K kg) gas constant rkice =2.04 ! W / (m K) ice conductivity hocond =2.5e6 ! J / kg heat of condensation/vap fusi =3.02e8 ! J / m^3 heat of fusion of ice fuss =1.10e8 ! J / m^3 heat of fusion of snow fice_max =.995 ! -- maximum fractional ice cover tice_m =273.05 ! K melting point of ice tsnw_m =273.15 ! K melting point of snow hice_min =.1 ! m minimum ice thickness epsmol =.622 ! -- molecular weight of h2o/dry air !albsa =12*.75 albsa =albs ! Seasonal albedo -- set to fixed value albs above sice = 6. ! per mil salinity of seaice rksnw =.31 ! W / (m K) snow conductivity cpair =1004. ! J / (K kg) specific heat of dry air cpsw =3987. ! J / (K kg) specific heat of seawater stefanb =5.67e-8 ! W / (m^2 K^4) stefan-boltzman constant aice =9.5 ! -- vapor pressure parameters bice =7.66 ! K .. awater =7.5 ! -- .. bwater =35.86 ! K .. t0deg =273.15 ! K zero deg celcius in K #if defined (ALBSNW_EVOL) albs_min =.71 ! -- minimum snow albedo albs_max =.85 ! -- maximum snow albedo taua =.008 ! -- constant for linear decrease of albedo for dry snow tauf =.24 ! -- constant for exponential decrease of wet snow due to wet metamorphism swe_newalb =0.002 ! m snowfall water equivalent depth necessary to refresh albedo back to maximum value hsnw_lim =.02 ! m limit snow depth to compute snow cover fraction (CICE formulation) #endif #if defined (SSNOWD_ICE) cv_snw =.68 ! -- coefficient of variation of the lognormal snow depth distribution ! It depends on the ice type and age. Based on SHEBA dataset for sea ice. #endif end subroutine icedat C C C subroutine iceinit !use mod_forcing_nersc implicit none integer :: i,j real :: flagi1,flagi2,tfrz, flagi, flaghemi include 'common_blocks.h' if (mnproc==1) then write(lp,*) 'Calling iceinit ' end if !$OMP PARALLEL DO PRIVATE(j,i,flagi,flaghemi,tfrz) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin do i=1-margin,ii+margin tfrz=273.216-.057*saln(i,j,1,1) tfrz=tfrz-t0deg c ! 1 if temperature < freezing point + .3 ! 0 otherwise flagi=(1.+sign(1.,tfrz+.3 -temp(i,j,1,1)))/2. c ! 0 in -50 < lat < 50 flagi=flagi*(1.+sign(1.,abs(plat(i,j))-50.))/2. c ! Hemisphere flag, 1 if plat > 0 flaghemi=(1.+sign(1.,plat(i,j)))/2. c ! Initial ice thickness - different for N/S hemispheres hicem(i,j)=flagi*(flaghemi*2.5 + (1.-flaghemi)*1.0) c ficem(i,j)=flagi*.95 hsnwm(i,j)=0. ticem(i,j)=temp(i,j,1,1)-flagi*4.+273.15 tsrfm(i,j)=ticem(i,j) qfrz(i,j)=0.0 pemnp(i,j)=0. #if defined (TEST_ICE_AGE) fy_age (i,j)=0. fy_frac (i,j)=0. rdg_frac(i,j)=0. #endif #if defined (ALBSNW_EVOL) albsnwm(i,j)=(albs_min+albs_max)*.5 #endif #if defined(SSNOWD_ICE) hmelt(i,j)=0. hprcp(i,j)=hsnwm(i,j) #endif enddo enddo !$OMP END PARALLEL DO c tauxice=0. tauyice=0. MIZ_MASK = 0.0!!test variable #if defined (ICE_DYN_DIAG) if ((itest.gt.0).and.(jtest.gt.0)) then !! on the tile that (itest,jest ) is on; !! - get mnproc for this tile !! - and the relative values of (itest,jtest) !! - sometimes common_blocks.h can't be used mnproc_2 = mnproc itest_2 = itest!!relative value of itest when in its own tile jtest_2 = jtest end if #endif /* ICE_DYN_DIAG */ end subroutine iceinit c subroutine clat_turb c c --- ------------------------------------------------------------------- c --- read numerical values of the latent transfer coefficient based c --- on the tabelled values in Isemer et al, J Clim., p. 1180, 1989 c --- n o t e : c --- i-index gives the wind at 10m height c --- from 0 to 30 m/s in intervals of 2 m/s c --- j-index gives the virtual air-sea temperature difference at 10m height c --- from -7 to +7 deg c in intervals of .5 deg c c --- for all but the equatorial and sub-equatorial waters, virtual temp is c --- close to real temp (see gill, 1982, p. 41) c --- ------------------------------------------------------------------- c use mod_xc implicit none REAL d_wind,d_temp,clat1(16,29) integer skip integer i,j logical :: ex c c dimension clat1(16,29) c c --- read data file inquire(exist=ex,file='iwh_tabulated.dat') if (.not.ex) then if (mnproc==1) & write(lp,*) 'clat_turb: iwh_tabulated.dat does not exist' call xcstop ('(clat_turb)') stop '(clat_turb)' end if open(19,file='iwh_tabulated.dat',STATUS= 'OLD') if (mnproc==1) print *,'reading wh_tabulated.dat' do 100 j=1,29 do 100 i=1,16 read(19,*) d_wind,d_temp,clat1(i,j) clat1(i,j)=clat1(i,j)*1.e-3 100 continue close(19) c c --- let i=1 represents wind speeds in the interval [0, 2) m/s c --- i=2 represents wind speeds in the interval [2, 4) m/s ... c do 111 j=1,29 do 110 i=1,15 110 clat(i,j)=(clat1(i,j)+clat1(i+1,j))*.5 111 clat(16,j)=clat(15,j) c c c --- let j=1 represents temp differences in the interval [-7, -6.5) deg c c --- j=2 represents temp differences in the interval [-6.5, -6) deg c... c do 121 i=1,16 do 120 j=1,28 120 clat(i,j)=(clat1(i,j)+clat1(i,j+1))*.5 121 clat(i,29)=clat(i,28) c end subroutine clat_turb c c c real function icevolume() implicit none include 'common_blocks.h' integer i,j,l real,dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: vol real*8 :: icevol8 vol=0.0 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP&SCHEDULE(STATIC,jblk) do j=1,jj do l=1,isp(j) do i=ifp(j,l),ilp(j,l) vol(i,j)=ficem(i,j)*hicem(i,j)*scpx(i,j)*scpy(i,j) end do end do end do !$OMP END PARALLEL DO !Gather tot ice volume from tiles call xcsum(icevol8,vol,ip) icevolume=icevol8 end function icevolume c c real function iceextent() implicit none include 'common_blocks.h' integer i,j,l real extent(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) real*8 :: iceext8 extent=0.0 !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP&SCHEDULE(STATIC,jblk) do 100 j=1,jj do 100 l=1,isp(j) do 100 i=ifp(j,l),ilp(j,l) extent(i,j)=ficem(i,j)*scpx(i,j)*scpy(i,j) 100 continue !$OMP END PARALLEL DO c --- Gather ice extent from tiles to get tot extent call xcsum(iceext8,extent,ip) iceextent=iceext8 end function iceextent c --- Ice volume/ice area diagnostics subroutine diag_icevol(dtime) use mod_year_info, only: daysinyear use mod_hycom_nersc #if defined(ICESTATE) use mod_icestate , only : icestate use mod_icestate_diag #endif implicit none real*8, intent(in) :: dtime integer :: iyear, iday, ihour, diy integer*8 :: ihour2,isec real :: icevol,iceext real, dimension(itdm,jtdm) :: modlat include 'common_blocks.h' c call forday(dtime, yrflag, iyear,iday,ihour) call xcaget(modlat,plat,0) ihour2=nint(dtime*86400.d0,kind=8)/3600 isec =nint(dtime*86400.d0,kind=8)-ihour2*3600 diy=daysinyear(iyear,yrflag) if (mod(ihour,6)==0 .and. isec < baclin) then #if defined (ICESTATE) icevol=icestate_volume() iceext=icestate_extent(.15) #else icevol=icevolume() iceext=iceextent() #endif if (mnproc==1) then open(10,file=rungen//'icevolume.dat',position='append') write(10,'(i5,g15.9,3g15.9)') & imem, iyear+(iday+(ihour+isec/3600.d0)/24.d0)/diy, & dtime,icevol*1.0d-9,iceext*1.0d-12 close(10) end if if (icevol < 1.0d12 .and. maxval(modlat) > 85) then call xcstop('(diag_icevol:Arctic Ice disappeared)') end if end if end subroutine diag_icevol #if defined (ICESTATE) ! =================================================================================== ! =========================== icestate2ice ========================================== ! =================================================================================== ! Routine transfers values from icestate to 'common_ice.h' variables. ! For diagnostic purposes (HYCOM prtsol) subroutine icestate2ice use mod_icestate , only:icestate implicit none integer i,j,l C$OMP PARALLEL DO PRIVATE(j,l,i) C$OMP&SCHEDULE(STATIC,1) do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) ficem(i,j) = sum(icestate(i,j)%ice%fice) hicem(i,j) = sum(icestate(i,j)%ice%fice*icestate(i,j)%ice%hice) hsnwm(i,j) = sum(icestate(i,j)%ice%fice*icestate(i,j)%ice%hsnw) ticem(i,j) = sum(icestate(i,j)%ice%fice*icestate(i,j)%ice%tsrf) ticem(i,j) = ticem(i,j) / (ficem(i,j) + 1e-6) hicem(i,j) = hicem(i,j) / (ficem(i,j) + 1e-6) hsnwm(i,j) = hsnwm(i,j) / (ficem(i,j) + 1e-6) tsrfm(i,j) = ticem(i,j)*ficem(i,j) + (1.-ficem(i,j))* & icestate(i,j)%tml end do end do end do C$OMP END PARALLEL DO end subroutine #endif c Ice advection routine c c c Ice advection routine c c subroutine iceadv(h,u,v,scuy,scvx,scp2i,scp2,dt) ! ! --- ------------------------------------------------------------------ ! --- Advection is done with flux limited 3rd order WENO in space and ! --- 2nd order Runge-Kutta in time ! --- ------------------------------------------------------------------ ! use mod_xc ! implicit none ! real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & & h,u,v,scuy,scvx,scp2i,scp2 real dt ! real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: sao,hp real dtm integer i,j,l sao(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy)=0 ! ! --- Use a modified time step since velocities are in m/s while scale ! --- factors are in cm ! dtm=dt*1.e2 dtm=dt ! --- Prediction step call weno3pd(h,sao,u,v,scuy,scvx,scp2i,scp2,dtm) !$OMP PARALLEL DO margin=nbdy do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) hp(i,j)=h(i,j)+dtm*sao(i,j) enddo enddo enddo !$OMP END PARALLEL DO ! --- Correction step call weno3pd(hp,sao,u,v,scuy,scvx,scp2i,scp2,dtm) !$OMP PARALLEL DO do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) h(i,j)=.5*(h(i,j)+hp(i,j)+dtm*sao(i,j)) enddo enddo enddo !$OMP END PARALLEL DO return end subroutine iceadv subroutine weno3pd(g,sao,u,v,scuy,scvx,scp2i,scp2,dt) ! ! --- ------------------------------------------------------------------ ! --- By a weighted essentially non-oscillatory scheme with up to 3th ! --- order accuracy, obtain the spatial advective operator of a ! --- 2-dimensional field defined at the scalar points of a C-grid. The ! --- fluxes are limited to make the scheme positive definite. ! --- Advective velocities in the i- and j-direction are defined at u- ! --- and v-points, respectively. ! --- ------------------------------------------------------------------ ! use mod_xc ! implicit none ! real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & & g,sao,u,v,scuy,scvx,scp2i,scp2 real dt ! real cq00,cq01,cq10,cq11,ca0,ca1,eps parameter (cq00=-1./2.,cq01= 3./2., & & cq10= 1./2.,cq11= 1./2., & & ca0=1./3.,ca1=2./3., & & eps=1.e-12) ! real, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & & ful,fuh,fvl,fvh,gt real q0,q1,a0,a1,q integer i,j,l,im1,im2,ip1,jm1,jm2,jp1 ! ! --- Compute grid cell boundary fluxes. Split in a low order flux ! --- (donor cell) and a high order correction flux. ! !$OMP PARALLEL DO do j=0,jj+2 do i=0,ii+2 ful(i,j)=0. fuh(i,j)=0. fvl(i,j)=0. fvh(i,j)=0. enddo enddo !$OMP END PARALLEL DO ! ! call xctilr(g, 1,1, 3,3, halo_ps) ! !$OMP PARALLEL DO PRIVATE(im1,im2,q0,q1,a0,a1,ip1) do j=0,jj+1 do l=1,isu(j) do i=max(0,ifu(j,l)),min(ii+2,ilu(j,l)) c do i=1,ii im1=i-1 ! if (u(i,j).gt.0.) then im2=im1-iu(im1,j) ! q0=cq00*g(im2,j)+cq01*g(im1,j) q1=cq10*g(im1,j)+cq11*g(i ,j) ! a0=ca0 a1=ca1*(abs(g(im2,j)-g(im1,j))+eps) & & /(abs(g(im1,j)-g(i ,j))+eps) ! ful(i,j)=u(i,j)*g(im1,j)*scuy(i,j) ! else ip1=i+iu(i+1,j) ! q0=cq11*g(im1,j)+cq10*g(i ,j) q1=cq01*g(i ,j)+cq00*g(ip1,j) ! a0=ca1 a1=ca0*(abs(g(im1,j)-g(i ,j))+eps) & & /(abs(g(i ,j)-g(ip1,j))+eps) ! ful(i,j)=u(i,j)*g(i ,j)*scuy(i,j) ! endif ! fuh(i,j)=u(i,j)*(a0*q0+a1*q1)/(a0+a1)*scuy(i,j)-ful(i,j) ! enddo enddo enddo !$OMP END PARALLEL DO ! !$OMP PARALLEL DO PRIVATE(jm1,q0,q1,a0,a1,jm2,jp1) do j=0,jj+2 jm1=j-1 do l=1,isv(j) do i=max(0,ifv(j,l)),min(ii+1,ilv(j,l)) c do i=1,ii ! if (v(i,j).gt.0.) then jm2=jm1-iv(i,jm1) ! q0=cq00*g(i,jm2)+cq01*g(i,jm1) q1=cq10*g(i,jm1)+cq11*g(i,j ) ! a0=ca0 a1=ca1*(abs(g(i,jm2)-g(i,jm1))+eps) & & /(abs(g(i,jm1)-g(i,j ))+eps) ! fvl(i,j)=v(i,j)*g(i,jm1)*scvx(i,j) ! else jp1=j+iv(i,j+1) ! q0=cq11*g(i,jm1)+cq10*g(i,j ) q1=cq01*g(i,j )+cq00*g(i,jp1) ! a0=ca1 a1=ca0*(abs(g(i,jm1)-g(i,j ))+eps) & & /(abs(g(i,j )-g(i,jp1))+eps) ! fvl(i,j)=v(i,j)*g(i,j )*scvx(i,j) ! endif ! fvh(i,j)=v(i,j)*(a0*q0+a1*q1)/(a0+a1)*scvx(i,j)-fvl(i,j) ! enddo enddo enddo !$OMP END PARALLEL DO ! ! --- Update field with low order fluxes. !$OMP PARALLEL DO do j=0,jj+1 do l=1,isp(j) do i=max(0,ifp(j,l)),min(ii+1,ilp(j,l)) ! do i=1,ii gt(i,j)=g(i,j)-dt*(ful(i+1,j)-ful(i,j) & & +fvl(i,j+1)-fvl(i,j))*scp2i(i,j) enddo enddo enddo !$OMP END PARALLEL DO ! ! --- Obtain fluxes with limited high order correction fluxes. q=.25/dt !$OMP PARALLEL DO do j=1,jj do l=1,isu(j) do i=max(1,ifu(j,l)),min(ii+1,ilu(j,l)) c do i=1,ii fuh(i,j)=ful(i,j)+max(-q*gt(i ,j)*scp2(i ,j), & & min( q*gt(i-1,j)*scp2(i-1,j),fuh(i,j))) enddo enddo enddo !$OMP END PARALLEL DO !$OMP PARALLEL DO do j=1,jj+1 do l=1,isv(j) do i=max(1,ifv(j,l)),min(ii,ilv(j,l)) c do i=1,ii fvh(i,j)=fvl(i,j)+max(-q*gt(i,j )*scp2(i,j ), & & min( q*gt(i,j-1)*scp2(i,j-1),fvh(i,j))) enddo enddo enddo !$OMP END PARALLEL DO ! ! --- Compute the spatial advective operator. !$OMP PARALLEL DO do j=1,jj do l=1,isp(j) do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) c do i=1,ii sao(i,j)=-(fuh(i+1,j)-fuh(i,j)+fvh(i,j+1)-fvh(i,j))*scp2i(i,j) enddo enddo enddo !$OMP END PARALLEL DO ! return end subroutine weno3pd #if defined(SSNOWD_ICE) subroutine fpsolver(cv,D_a,D_ave,D_m) use m_get_erfc c This uses a "fixed-point" iteration algorithm to solve for the c depth of melt when the cell-average snow depth is known. implicit none integer maxiter,i real cv,D_a,D_m,sca_Dm,zeta,D_ave,xlambda,z_Dm,tol,D_m_old c Define the melt-depth tolerance to be 0.1 mm. tol = 1.0e-4 maxiter = 10 c Set the initial guess to a small number (this does not seem to c affect the iterations required to gain convergence, and c starting with a large number can lead to divergence of the c solution). D_m_old = 1.0e-9 zeta = sqrt(log(1.0 + max(cv,0.001)**2)) xlambda = log(D_a) - 0.5 * zeta**2 do i=1,maxiter z_Dm = (log(D_m_old) - xlambda) / zeta sca_Dm = 0.5 * get_erfc(z_Dm/sqrt(2.0)) D_m = (0.5 * exp(xlambda + 0.5*zeta**2) * & get_erfc((z_Dm - zeta)/sqrt(2.0)) - D_ave) / sca_Dm c print *, i,cv,D_m,sca_Dm,D_a,D_ave if (abs(D_m - D_m_old).lt.tol) return D_m_old = D_m return end subroutine fpsolver #endif/* SSNOWD_ICE */ #if defined (ICE_DYN_DIAG) && defined (ICE_DYN_DIAG_DUMP) subroutine ice_dyn_diag_dump(dtime) !!routine to dump the following variables periodically: !! iceu !! icev !! strainI !! strainII !! stressp !! stressm !! stress12 !! pice !! dfloe !! MIZ_mask use mod_xc use mod_hycom_nersc!->ncdraft use netcdf #if defined (WAVES) use mod_common_wavesice, only: dfloe #endif implicit none real, intent(in) :: dtime !!for determining when to dump ice stuff real*8 ,save :: dtime_icedump_last!!this is set to exactly zero on first call logical :: DO_ICEDUMP real*8 :: icedump_timestep !!time between dumps (in days) & ,dtime_icedump & ,diff_dtime_icedump !!for netcdf character(len=80) :: ncfil integer :: iyear, ihour, iday character(len=4) :: cyy character(len=3) :: cdd,cdfrac character(len=12) :: dtime_human!!yyyy_ddd.xxx, where dd starts at 0 character(len=*), parameter :: outdir = './ICE_DYN_DIAG/' !! real, dimension(itdm,jtdm) :: fld include 'common_blocks.h' icedump_timestep = icedump_freq/(24.0)!!dump frequency in days if (DO_INIT_ICEDUMP) then DO_INIT_ICEDUMP = .false. DO_ICEDUMP = .false.!!don't dump at time 0 if (dump_time0==.false.) then dtime_icedump_last = dtime else DO_ICEDUMP = .true.!!do dump at time 0 dtime_icedump = dtime !!do dump next timestep since some fields aren't in restart file !dtime_icedump_last = dtime-icedump_timestep end if if (mnproc==1) then print*,'TWice_diag: initialising ice_dyn_diag_dump' !print*,'TWice_diag',dtime_icedump_last end if else dtime_icedump = floor(dtime/icedump_timestep) & *icedump_timestep !! diff_dtime_icedump = (dtime_icedump - dtime_icedump_last) DO_ICEDUMP = (diff_dtime_icedump.ge.icedump_timestep) end if if (DO_ICEDUMP) then dtime_icedump_last = dtime_icedump !!set name of nc file call forday(dtime_icedump,yrflag, iyear,iday,ihour) write(cyy ,'(i4.4)') iyear write(cdd,'(i3.3)') iday-1!!start at iday=0 write(cdfrac,'(i3.3)') int(1.0e3*(ihour/24.0)) dtime_human = cyy//'_'//cdd//'.'//cdfrac ncfil = outdir//'/ice_dyn_diag_'//dtime_human//'.nc' if (mnproc==1) then print*,'TWice_diag: making ncfil: ',ncfil,dtime,dtime_human print*,'TWice_diag: time: ',dtime,dtime_human end if !!dump lon/lat call xcaget(fld,plon,0) call ncdraft(trim(ncfil),fld,'lon','clobber') call xcaget(fld,plat,0) call ncdraft(trim(ncfil),fld,'lat','') !!dump icef/iceh call xcaget(fld,ficem,0) call ncdraft(trim(ncfil),fld,'ficem','') call xcaget(fld,hicem,0) call ncdraft(trim(ncfil),fld,'hicem','') !!dump iceu/icev call xcaget(fld,iceu,0) call ncdraft(trim(ncfil),fld,'iceu','') call xcaget(fld,icev,0) call ncdraft(trim(ncfil),fld,'icev','') !!dump strainI/strainII call xcaget(fld,strainI,0) call ncdraft(trim(ncfil),fld,'strI','') call xcaget(fld,strainII,0) call ncdraft(trim(ncfil),fld,'strII','') !!dump stressp/stressm/stress12 call xcaget(fld,stressp,0) call ncdraft(trim(ncfil),fld,'sigp','') call xcaget(fld,stressm,0) call ncdraft(trim(ncfil),fld,'sigm','') call xcaget(fld,stress12,0) call ncdraft(trim(ncfil),fld,'sig12','') call xcaget(fld,pice,0) call ncdraft(trim(ncfil),fld,'pice','') !!dump MIZ_MASK call xcaget(fld,MIZ_MASK,0) call ncdraft(trim(ncfil),fld,'MIZ_MASK','') #if defined (WAVES) !!dump dfloe call xcaget(fld,dfloe,0) call ncdraft(trim(ncfil),fld,'dfloe','') #endif if (mnproc==1) then print*,'TWice_diag: finished making ncfil: ',ncfil end if else if (mnproc==1) then ncfil = 'ice_dyn_diag_' print*,'TWid2: not making ncfil; time: ',dtime,dtime_human end if end if end subroutine ice_dyn_diag_dump #endif/*ICE_DYN_DIAG && ICE_DYN_DIAG_DUMP*/ end module mod_common_ice