      subroutine improvedclouds(ngrid,nlay,ptimestep,
     &             pplay,pt,pdt,
     &             pq,pdq,pdqcloud,pdtcloud,
     &             nq,tauscaling)
! to use  'getin'
      USE ioipsl_getincom
      USE updaterad
      implicit none
      
      
c------------------------------------------------------------------
c  This routine is used to form clouds when a parcel of the GCM is
c    saturated. It includes the ability to have supersaturation, a
c    computation of the nucleation rates, growthrates and the
c    scavenging of dust particles by clouds.
c  It is worth noting that the amount of dust is computed using the
c    dust optical depth computed in aeropacity.F. That's why
c    the variable called "tauscaling" is used to convert
c    pq(dust_mass) and pq(dust_number), which are relative
c    quantities, to absolute and realistic quantities stored in zq.
c    This has to be done to convert the inputs into absolute
c    values, but also to convert the outputs back into relative
c    values which are then used by the sedimentation and advection
c    schemes.

c  Authors: J.-B. Madeleine, based on the work by Franck Montmessin
c           (October 2011)
c           T. Navarro, debug,correction, new scheme (October-April 2011)
c           A. Spiga, optimization (February 2012)
c------------------------------------------------------------------
#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "callkeys.h"
#include "tracer.h"
#include "comgeomfi.h"
#include "dimradmars.h"
#include "microphys.h"
#include "conc.h"
c------------------------------------------------------------------
c     Inputs:

      INTEGER ngrid,nlay
      integer nq                 ! nombre de traceurs
      REAL ptimestep             ! pas de temps physique (s)
      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
            
      REAL pt(ngrid,nlay)        ! temperature at the middle of the
                                 !   layers (K)
      REAL pdt(ngrid,nlay)       ! tendance temperature des autres
                                 !   param.
      REAL pq(ngrid,nlay,nq)     ! traceur (kg/kg)
      REAL pdq(ngrid,nlay,nq)    ! tendance avant condensation
                                 !   (kg/kg.s-1)
      REAL tauscaling(ngridmx)     ! Convertion factor for qdust and Ndust

c     Outputs:
      REAL pdqcloud(ngrid,nlay,nq) ! tendance de la condensation
                                   !   H2O(kg/kg.s-1)
      REAL pdtcloud(ngrid,nlay)    ! tendance temperature due
                                   !   a la chaleur latente

c------------------------------------------------------------------
c     Local variables:

      LOGICAL firstcall
      DATA firstcall/.true./
      SAVE firstcall

      REAL*8   derf ! Error function
      !external derf
     
      INTEGER ig,l,i
      
      REAL zq(ngridmx,nlayermx,nqmx)  ! local value of tracers
      REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers
      REAL zt(ngridmx,nlayermx)       ! local value of temperature
      REAL zqsat(ngridmx,nlayermx)    ! saturation
      REAL lw                         !Latent heat of sublimation (J.kg-1) 
      REAL cste
      REAL dMice           ! mass of condensed ice
!      REAL sumcheck
      REAL*8 ph2o          ! Water vapor partial pressure (Pa)
      REAL*8 satu          ! Water vapor saturation ratio over ice
      REAL*8 Mo,No
      REAL*8 Rn, Rm, dev2, n_derf, m_derf
      REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
      REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin

      REAL*8 sig      ! Water-ice/air surface tension  (N.m)
      EXTERNAL sig

      REAL dN,dM
      REAL rate(nbin_cld)  ! nucleation rate
      REAL seq

      REAL rice(ngrid,nlay)      ! Ice mass mean radius (m)
                                 ! (r_c in montmessin_2004)
      REAL rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
      REAL rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)

      REAL res      ! Resistance growth
      

c     Parameters of the size discretization
c       used by the microphysical scheme
      DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
      DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
      DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6
                                           ! Minimum boundary radius (m)
      DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
      DOUBLE PRECISION vrat_cld ! Volume ratio
      DOUBLE PRECISION rb_cld(nbin_cld+1)! boundary values of each rad_cld bin (m)
      SAVE rb_cld
      DOUBLE PRECISION dr_cld(nbin_cld)   ! width of each rad_cld bin (m)
      DOUBLE PRECISION vol_cld(nbin_cld)  ! particle volume for each bin (m3)


      REAL sigma_ice ! Variance of the ice and CCN distributions
      SAVE sigma_ice


      
c----------------------------------      
c TESTS

      INTEGER countcells
      
      LOGICAL test_flag    ! flag for test/debuging outputs
      SAVE    test_flag    


      REAL satubf(ngrid,nlay),satuaf(ngrid,nlay) 
      REAL res_out(ngrid,nlay)
 

c------------------------------------------------------------------

      IF (firstcall) THEN
!=============================================================
! 0. Definition of the size grid
!=============================================================
c       rad_cld is the primary radius grid used for microphysics computation.
c       The grid spacing is computed assuming a constant volume ratio
c       between two consecutive bins; i.e. vrat_cld.
c       vrat_cld is determined from the boundary values of the size grid: 
c       rmin_cld and rmax_cld.
c       The rb_cld array contains the boundary values of each rad_cld bin.
c       dr_cld is the width of each rad_cld bin.

c       Volume ratio between two adjacent bins
   !     vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
   !     vrat_cld = exp(vrat_cld)
        vrat_cld = dlog(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
        vrat_cld = dexp(vrat_cld)
        write(*,*) "vrat_cld", vrat_cld

        rb_cld(1)  = rbmin_cld
        rad_cld(1) = rmin_cld
        vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld
   !     vol_cld(1) = 4./3. * pi * rmin_cld*rmin_cld*rmin_cld

        do i=1,nbin_cld-1
          rad_cld(i+1)  = rad_cld(i) * vrat_cld**(1./3.)
          vol_cld(i+1)  = vol_cld(i) * vrat_cld
        enddo
        
        do i=1,nbin_cld
          rb_cld(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) *
     &      rad_cld(i)
          dr_cld(i)  = rb_cld(i+1) - rb_cld(i)
        enddo
        rb_cld(nbin_cld+1) = rbmax_cld
        dr_cld(nbin_cld)   = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)

        print*, ' '
        print*,'Microphysics: size bin information:'
        print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
        print*,'-----------------------------------'
        do i=1,nbin_cld
          write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i),
     &      dr_cld(i)
        enddo
        write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
        print*,'-----------------------------------'

        do i=1,nbin_cld+1
!           rb_cld(i) = log(rb_cld(i))  
            rb_cld(i) = dlog(rb_cld(i))  !! we save that so that it is not computed
                                         !! at each timestep and gridpoint
        enddo

c       Contact parameter of water ice on dust ( m=cos(theta) )
c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!       mteta  = 0.95
        write(*,*) 'water_param contact parameter:', mteta

c       Volume of a water molecule (m3)
        vo1 = mh2o / dble(rho_ice)
c       Variance of the ice and CCN distributions
        sigma_ice = sqrt(log(1.+nuice_sed))
       
 
        write(*,*) 'Variance of ice & CCN distribs :', sigma_ice
        write(*,*) 'nuice for sedimentation:', nuice_sed
        write(*,*) 'Volume of a water molecule:', vo1


        test_flag = .false.
 
        firstcall=.false.
      END IF


!=============================================================
! 1. Initialisation
!=============================================================
      cste = 4*pi*rho_ice*ptimestep

      res_out(:,:) = 0
      rice(:,:) = 1.e-8

c     Initialize the tendencies
      pdqcloud(1:ngrid,1:nlay,1:nq)=0
      pdtcloud(1:ngrid,1:nlay)=0
    
c     Initialize the tendencies
      pdqcloud(1:ngrid,1:nlay,1:nq)=0
      pdtcloud(1:ngrid,1:nlay)=0
      
c     Initialize the tendencies
      pdqcloud(1:ngrid,1:nlay,1:nq)=0
      pdtcloud(1:ngrid,1:nlay)=0
      
      
      zt(1:ngrid,1:nlay) = 
     &      pt(1:ngrid,1:nlay) + 
     &      pdt(1:ngrid,1:nlay) * ptimestep

      zq(1:ngrid,1:nlay,1:nq) = 
     &      pq(1:ngrid,1:nlay,1:nq) + 
     &      pdq(1:ngrid,1:nlay,1:nq) * ptimestep
      
      
      WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
     &       zq(1:ngrid,1:nlay,1:nq) = 1.e-30

      zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
      
!=============================================================
! 2. Compute saturation
!=============================================================


      dev2 = 1. / ( sqrt(2.) * sigma_ice )

      call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
            
      countcells = 0

c     Main loop over the GCM's grid
      DO l=1,nlay
        DO ig=1,ngrid

c       Get the partial pressure of water vapor and its saturation ratio
        ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
        satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)

!=============================================================
! 3. Nucleation
!=============================================================

       IF ( satu .ge. 1. ) THEN         ! if there is condensation

        call updaterccn(zq(ig,l,igcm_dust_mass),
     &          zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))


c       Expand the dust moments into a binned distribution
        Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)   + 1.e-30
        No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
        Rn = rdust(ig,l)
        Rn = -dlog(Rn) 
        Rm = Rn - 3. * sigma_ice*sigma_ice  
        n_derf = derf( (rb_cld(1)+Rn) *dev2)
        m_derf = derf( (rb_cld(1)+Rm) *dev2)
        do i = 1, nbin_cld
          n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
          m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
          n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
          m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
          n_aer(i) = n_aer(i) + 0.5 * No * n_derf
          m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
        enddo
        
!        sumcheck = 0
!        do i = 1, nbin_cld
!          sumcheck = sumcheck + n_aer(i)
!        enddo
!        sumcheck = abs(sumcheck/No - 1)
!        if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
!          print*, "WARNING, No sumcheck PROBLEM"
!          print*, "sumcheck, No",sumcheck, No
!          print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
!          print*, "Dust binned distribution", n_aer
!        endif
!        
!        sumcheck = 0
!        do i = 1, nbin_cld
!          sumcheck = sumcheck + m_aer(i)
!        enddo
!        sumcheck = abs(sumcheck/Mo - 1)
!        if ((sumcheck .gt. 1e-5) .and.  (1./Rn .gt. rmin_cld)) then
!          print*, "WARNING, Mo sumcheck PROBLEM"
!          print*, "sumcheck, Mo",sumcheck, Mo
!          print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
!          print*, "Dust binned distribution", m_aer
!        endif

  
c       Get the rates of nucleation
        call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
        
        dN = 0.
        dM = 0.
        do i = 1, nbin_cld
          n_aer(i) = n_aer(i)/( 1. + rate(i)*ptimestep)
          m_aer(i) = m_aer(i)/( 1. + rate(i)*ptimestep)
          dN       = dN + n_aer(i) * rate(i) * ptimestep
          dM       = dM + m_aer(i) * rate(i) * ptimestep
        enddo


c       Update Dust particles
        zq(ig,l,igcm_dust_mass)   = 
     &  zq(ig,l,igcm_dust_mass)   - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10) 
        zq(ig,l,igcm_dust_number) = 
     &  zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
c       Update CCNs
        zq(ig,l,igcm_ccn_mass)   = 
     &  zq(ig,l,igcm_ccn_mass)   + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
        zq(ig,l,igcm_ccn_number) = 
     &  zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
        
        ENDIF ! of is satu >1

!=============================================================
! 4. Ice growth: scheme for radius evolution
!=============================================================

c We trigger crystal growth if and only if there is at least one nuclei (N>1).
c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.

       IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth 

     
        call updaterice_micro(zq(ig,l,igcm_h2o_ice),
     &    zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
     &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))

        No   = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30

c       saturation at equilibrium
c       rice should not be too small, otherwise seq value is not valid
        seq  = exp(2.*sig(zt(ig,l))*mh2o / (rho_ice*rgp*zt(ig,l)*
     &             max(rice(ig,l),1.e-7)))
        
c       get resistance growth
        call growthrate(zt(ig,l),pplay(ig,l),
     &             real(ph2o/satu),rice(ig,l),res)

        res_out(ig,l) = res

ccccccc  implicit scheme of mass growth

        dMice =
     & (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))
     &   /(res*zqsat(ig,l)/(cste*No*rice(ig,l)) + 1.)


! With the above scheme, dMice cannot be bigger than vapor, 
! but can be bigger than all available ice.
       dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
       dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...

       zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
       zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice


       countcells = countcells + 1 
       
       ! latent heat release
       lw=(2834.3-0.28*(zt(ig,l)-To)-
     &     0.004*(zt(ig,l)-To)*(zt(ig,l)-To))*1.e+3
       pdtcloud(ig,l)= dMice*lw/cpp/ptimestep
          
          
     
!=============================================================
! 5. Dust cores released, tendancies, latent heat, etc ...
!=============================================================

c         If all the ice particles sublimate, all the condensation
c         nuclei are released:
          if (zq(ig,l,igcm_h2o_ice).le.1.e-8) then
          
c           Water
            zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap) 
     &                            + zq(ig,l,igcm_h2o_ice)
            zq(ig,l,igcm_h2o_ice) = 0.
c           Dust particles
            zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
     &                              + zq(ig,l,igcm_ccn_mass)
            zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
     &                                + zq(ig,l,igcm_ccn_number)
c           CCNs
            zq(ig,l,igcm_ccn_mass) = 0.
            zq(ig,l,igcm_ccn_number) = 0.

          endif
         
          ENDIF !of if Nccn>1
         
        ENDDO ! of ig loop
      ENDDO ! of nlayer loop
      
      
      ! Get cloud tendencies
        pdqcloud(1:ngrid,1:nlay,igcm_h2o_vap) =
     &   (zq(1:ngrid,1:nlay,igcm_h2o_vap) - 
     &    zq0(1:ngrid,1:nlay,igcm_h2o_vap))/ptimestep
        pdqcloud(1:ngrid,1:nlay,igcm_h2o_ice) =
     &   (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
     &    zq0(1:ngrid,1:nlay,igcm_h2o_ice))/ptimestep
        pdqcloud(1:ngrid,1:nlay,igcm_ccn_mass) =
     &   (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
     &    zq0(1:ngrid,1:nlay,igcm_ccn_mass))/ptimestep
        pdqcloud(1:ngrid,1:nlay,igcm_ccn_number) =
     &   (zq(1:ngrid,1:nlay,igcm_ccn_number) -
     &    zq0(1:ngrid,1:nlay,igcm_ccn_number))/ptimestep
     
      if (scavenging) then
      
        pdqcloud(1:ngrid,1:nlay,igcm_dust_mass) =
     &   (zq(1:ngrid,1:nlay,igcm_dust_mass) -
     &    zq0(1:ngrid,1:nlay,igcm_dust_mass))/ptimestep
        pdqcloud(1:ngrid,1:nlay,igcm_dust_number) =
     &   (zq(1:ngrid,1:nlay,igcm_dust_number) -
     &    zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep
          
      endif
     
     
     
!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 
!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 
!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
      IF (test_flag) then
      
!       error2d(:) = 0.
       DO l=1,nlay
       DO ig=1,ngrid
!         error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
          satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l) 
          satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l) 
       ENDDO
       ENDDO

      print*, 'count is ',countcells, ' i.e. ',
     &     countcells*100/(nlay*ngrid), '% for microphys computation'

!      IF (ngrid.ne.1) THEN ! 3D
!         call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
!     &                    satu_out)
!         call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
!     &                    dM_out)
!         call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
!     &                    dN_out)
!         call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
!     &                    error2d)
!         call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
!     &                    zqsat)
!      ENDIF

!      IF (ngrid.eq.1) THEN ! 1D
!         call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
!     &                    error_out)
         call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1,
     &                    res_out)
         call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1,
     &                    satubf)
         call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1,
     &                    satuaf)
         call WRITEdiagfi(ngrid,"vapbf","h2ovap before","kg/kg",1,
     &                    zq0(1,:,igcm_h2o_vap))
         call WRITEdiagfi(ngrid,"vapaf","h2ovap after","kg/kg",1,
     &                    zq(1,:,igcm_h2o_vap))
         call WRITEdiagfi(ngrid,"icebf","h2oice before","kg/kg",1,
     &                    zq0(1,:,igcm_h2o_ice))
         call WRITEdiagfi(ngrid,"iceaf","h2oice after","kg/kg",1,
     &                    zq(1,:,igcm_h2o_ice))
         call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1,
     &                    zq0(1,:,igcm_ccn_number))
         call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1,
     &                    zq(1,:,igcm_ccn_number))
c         call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
c     &                    gr_out)
c         call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
c     &                    rate_out)
c         call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
c     &                    dM_out)
c         call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
c     &                    dN_out)
         call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1,
     &                    zqsat)
!         call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
!     &                    satu_out)
         call WRITEdiagfi(ngrid,"rice","ice radius","m",1,
     &                    rice)
!         call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
!     &                    rdust)
!         call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
!     &                    rsedcloud)
!         call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
!     &                    rhocloud)
!      ENDIF
      
      ENDIF ! endif test_flag
!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 
!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 
!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 

      return
      end
      
      
      
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
c It is an analytical solution to the ice radius growth equation, 
c with the approximation of a constant 'reduced' cunningham correction factor 
c (lambda in growthrate.F) taken at radius req instead of rice    
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c      subroutine phi(rice,req,coeff1,coeff2,time)
c      
c      implicit none
c      
c      ! inputs
c      real rice ! ice radius
c      real req  ! ice radius at equilibirum
c      real coeff1  ! coeff for the log
c      real coeff2  ! coeff for the arctan
c
c      ! output      
c      real time
c      
c      !local
c      real var
c      
c      ! 1.73205 is sqrt(3)
c      
c      var = max(
c     &  abs(rice-req) / sqrt(rice*rice + rice*req  + req*req),1e-30)
c            
c       time = 
c     &   coeff1 * 
c     &   log( var )
c     & + coeff2 * 1.73205 *
c     &   atan( (2*rice+req) / (1.73205*req) )
c      
c      return
c      end
      
      
      
