source: trunk/LMDZ.MARS/libf/phymars/surfini.F @ 1711

      SUBROUTINE surfini(ngrid,piceco2,qsurf,psolaralb)
   ! to use  'getin'
      USE ioipsl_getincom
      use netcdf
      use tracer_mod, only: nqmx, noms
      use geometry_mod, only: longitude, latitude ! in radians
      use surfdat_h, only: watercaptag, frost_albedo_threshold,
     &                     albedo_h2o_ice, inert_h2o_ice, albedodat,
     &                     albedice, dryness
#ifndef MESOSCALE
      use mod_grid_phy_lmdz, only : klon_glo ! # of physics point on full grid
      use mod_phys_lmdz_para, only : is_master, gather, scatter
#endif
      USE comcstfi_h
      use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
      IMPLICIT NONE
c=======================================================================
c
c   creation des calottes pour l'etat initial
c
c=======================================================================
c-----------------------------------------------------------------------
c   Declarations:
c   -------------
#include "callkeys.h"
#include "datafile.h"

      integer,intent(in) :: ngrid ! number of atmospheric columns
      real,intent(in) :: piceco2(ngrid) ! CO2 ice thickness
      real,intent(inout) :: qsurf(ngrid,nqmx) ! tracer on surface (kg/m2)
      real,intent(out) :: psolaralb(ngrid,2) ! albedo

      INTEGER ig,icap,iq,alternate
      REAL icedryness ! ice dryness
      
      ! longwatercaptag is watercaptag. Trick for some compilers
      LOGICAL, DIMENSION(100000) :: longwatercaptag
      
! There are 3 different modes for ice distribution:
! icelocationmode = 1 ---> based on data from surface.nc
! icelocationmode = 2 ---> directly predefined for GCM resolutions 32x24 or 64x48
! icelocationmode = 3 ---> based on logical relations for latitude and longitude
! For visualisation : > /u/tnalmd/bin/watercaps gcm_txt_output_file
      INTEGER,SAVE :: icelocationmode = 2
       
       
      !in case icelocationmode == 1
      INTEGER i,j
      INTEGER     imd,jmd
      PARAMETER   (imd=360,jmd=180)
      REAL        zdata(imd,jmd)
      REAL        zelat,zelon 

#ifndef MESOSCALE
      INTEGER nb_ice(klon_glo,2)   ! number of counts | detected ice for GCM grid
#endif
      INTEGER latice(nbp_lat-1,2),lonice (nbp_lon,2) ! number of counts | detected ice along lat & lon axis

      REAL step,count,ratiolat

      INTEGER   ierr,nid,nvarid
      
      REAL,SAVE :: min_icevalue = 500.
      character(len=50) :: string = 'thermal'
      
      character (len=100) :: zedatafile

#ifdef MESOSCALE

      do ig=1,ngrid

         !write(*,*) "all qsurf to zero. dirty."
         do iq=1,nqmx
         qsurf(ig,iq)=0.  !! on jette les inputs GCM
                          !! on regle juste watercaptag
                          !! il faudrait garder les inputs GCM 
                          !! si elles sont consequentes
         enddo
         if ( ( latitude(ig)*180./pi .gt. 70. ) .and.
     .        ( albedodat(ig) .ge. 0.26   ) )  then
                 write(*,*)"outlier ",ig
                 watercaptag(ig)  = .true.
                 dryness(ig)      = 1.
                 albedodat(ig)    = albedo_h2o_ice  !! pour output 
         else
                 watercaptag(ig)  = .false.
                 dryness(ig)      = 1.
         endif

      enddo
#endif
! problem with nested precompiling flags

#ifndef MESOSCALE
      ! to handle parallel cases
#if CPP_PARA
      logical watercaptag_glo(klon_glo)
      real dryness_glo(klon_glo)
      real lati_glo(klon_glo)
      real long_glo(klon_glo)
#else
      logical watercaptag_glo(ngrid)
      real dryness_glo(ngrid)
      real lati_glo(ngrid)
      real long_glo(ngrid)
#endif
#endif

#ifndef MESOSCALE

c
c=======================================================================
! Initialize watercaptag (default is false)
      watercaptag_glo(:)=.false.

c     water ice outliers
c     ------------------------------------------

      IF ((water) .and. (caps)) THEN
     
c Perennial H20 north cap defined by watercaptag=true (allows surface to be
c hollowed by sublimation in vdifc).

c We might not want albedodat to be modified because it is used to write 
c restart files. Instead, albedo is directly modified when needed (i.e. 
c if we have watercaptag and no co2 ice), below and in albedocaps.F90

c       "Dryness coefficient" controlling the evaporation and
c        sublimation from the ground water ice (close to 1)
c        HERE, the goal is to correct for the fact
c        that the simulated permanent water ice polar caps
c        is larger than the actual cap and the atmospheric
c        opacity not always realistic.

         alternate = 0
         
         if (ngrid .ne. 1) then
           watercaptag(:) = .false.
           longwatercaptag(:) = .false.
         endif
         
         write(*,*) "surfini: Ice dryness ?"
         icedryness=1. ! default value
         call getin("icedryness",icedryness)
         write(*,*) "surfini: icedryness = ",icedryness
         dryness (:) = icedryness
         
      ! To be able to run in parallel, we work on the full grid
      ! and dispatch results afterwards

      ! start by geting latitudes and logitudes on full grid
      ! (in serial mode, this is just a copy)
      call gather(latitude,lati_glo)
      call gather(longitude,long_glo)

      if (is_master) then

        IF (ngrid .eq. 1) THEN ! special case for 1d --> do nothing
      
         print*, 'ngrid = 1, do no put ice caps in surfini.F'

        ELSE IF (icelocationmode .eq. 1) THEN
      
         print*,'Surfini: ice caps defined from surface.nc'
            
! This method detects ice as gridded value above min_icevalue in the field "string" from surface.nc
! Typically, it is for thermal inertia above 500 tiu.
! Two conditions are verified:
! 1. GCM ice caps are defined such as area is conserved for a given latitude
! (the approximation is that all points within the GCM latitude resolution have the same area).
! 2. caps are placed to fill the GCM points with the most detected ice first.
      

           
         zedatafile = trim(datafile)
 
        
         ierr=nf90_open(trim(zedatafile)//'/surface.nc',
     &   NF90_NOWRITE,nid)
     
         IF (ierr.NE.nf90_noerr) THEN
       write(*,*)'Error : cannot open file surface.nc '
       write(*,*)'(in phymars/surfini.F)'
       write(*,*)'It should be in :',trim(zedatafile),'/'
       write(*,*)'1) You can set this path in the callphys.def file:'
       write(*,*)'   datadir=/path/to/the/datafiles'
       write(*,*)'2) If necessary, surface.nc (and other datafiles)'
       write(*,*)'   can be obtained online on:'
       write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
       CALL ABORT
         ENDIF
      
      
         ierr=nf90_inq_varid(nid, string, nvarid)
         if (ierr.ne.nf90_noerr) then
          write(*,*) 'surfini error, cannot find ',trim(string)
          write(*,*) ' in file ',trim(zedatafile),'/surface.nc'
          write(*,*)trim(nf90_strerror(ierr))
          stop
         endif

         ierr=nf90_get_var(nid, nvarid, zdata)

         if (ierr.ne.nf90_noerr) then
          write(*,*) 'surfini: error failed loading ',trim(string)
          write(*,*)trim(nf90_strerror(ierr))
          stop
         endif
 
                     
         ierr=nf90_close(nid)
 

         nb_ice(:,1) = 1 ! default: there is no ice
         latice(:,1) = 1
         lonice(:,1) = 1
         nb_ice(:,2) = 0
         latice(:,2) = 0
         lonice(:,2) = 0
         !print*,'jjm,iim',jjm,iim ! jjm =  nb lati , iim = nb longi

         ! loop over the GCM grid - except for poles (ig=1 and ngrid)
         do ig=2,klon_glo-1
      
          ! loop over the surface file grid      
          do i=1,imd
           do j=1,jmd
             zelon = i - 180.
             zelat = 90. - j 
            if ((abs(lati_glo(ig)*180./pi-zelat).le.
     &           90./real(nbp_lat-1)) .and.
     &          (abs(long_glo(ig)*180./pi-zelon).le.
     &           180./real(nbp_lon))) then
              ! count all points in that GCM grid point
              nb_ice(ig,1) = nb_ice(ig,1) + 1
              if (zdata(i,j) > min_icevalue)
                 ! count all detected points in that GCM grid point
     &           nb_ice(ig,2) = nb_ice(ig,2) + 1
             endif
           enddo
          enddo  

        ! projection of nb_ice on GCM lat and lon axes
          latice(1+(ig-2)/nbp_lon,:) =
     &     latice(1+(ig-2)/nbp_lon,:) + nb_ice(ig,:)
          lonice(1+mod(ig-2,nbp_lon),:) = 
     &     lonice(1+mod(ig-2,nbp_lon),:) + nb_ice(ig,:) ! lonice is USELESS ...

         enddo ! of do ig=2,klon_glo-1
     

     
         ! special case for poles
         nb_ice(1,2)   = 1  ! ice prescribed on north pole
         latice(1,:)   = nb_ice(1,:)
         lonice(1,:)   = nb_ice(1,:)
         latice(nbp_lat-1,:) = nb_ice(ngrid,:)
         lonice(nbp_lon,:) = nb_ice(ngrid,:)
      
     
!      print*, 'latice TOT', latice(:,1)
!      print*, 'latice FOUND', latice(:,2)
!      print*, 'lonice TOT', lonice(:,1)
!      print*, 'lonice FOUND', lonice(:,2)
      
!      print*, 'lat ratio', int(real(latice(:,2))/real(latice(:,1))*iim)
!      print*, 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))*jjm)
      
!      print*,''
!      print*,'sum lat', sum(latice(:,1)), sum(lonice(:,1))
!      print*,'sum lon', sum(latice(:,2)), sum(lonice(:,2))
      
    
         ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
         do i=1,(nbp_lat-1)/2
          step  = 1. ! threshold to add ice cap
          count = 0. ! number of ice GCM caps at this latitude
          ! ratiolat is the ratio of area covered by ice within this GCM latitude range
          ratiolat  = real(latice(i,2))/real(latice(i,1))
          !print*,'i',i,(i-1)*iim+2,i*iim+1
     
          ! put ice caps while there is not enough ice,
          ! as long as the threshold is above 20%
          do while ((count.le.ratiolat*nbp_lon).and.(step.ge.0.2))
           count = 0.
           ! loop over GCM longitudes
           do j=1,nbp_lon
            ! if the detected ice ratio in the GCM grid point 
            ! is more than 'step', then add ice
            if (real(nb_ice((i-1)*nbp_lon+1+j,2)) 
     &        / real(nb_ice((i-1)*nbp_lon+1+j,1)) .ge. step) then
                  watercaptag_glo((i-1)*nbp_lon+1+j) = .true.
                  count = count + 1
            endif
           enddo ! of do j=1,nbp_lon
           !print*, 'step',step,count,ratiolat*nbp_lon
           step = step - 0.01
          enddo ! of do while
          !print*, 'step',step,count,ratiolat*nbp_lon

         enddo ! of do i=1,jjm/2
            

        ELSE IF (icelocationmode .eq. 2) THEN
      
         print*,'Surfini: predefined ice caps'
      
         if ((nbp_lon.eq.32).and.((nbp_lat-1).eq.24)) then ! 32x24
           
          print*,'water ice caps distribution for 32x24 resolution'
          longwatercaptag(1:9)    = .true. ! central cap - core
          longwatercaptag(26:33)  = .true. ! central cap
          longwatercaptag(1:33)  = .true. ! central cap
          longwatercaptag(56)  = .true. ! central cap
          longwatercaptag(58)  = .true. ! central cap
          longwatercaptag(60)  = .true. ! central cap
          longwatercaptag(62)  = .true. ! central cap
          longwatercaptag(64)  = .true. ! central cap
!---------------------   OUTLIERS  ----------------------------

         else if ((nbp_lon.eq.64).and.((nbp_lat-1).eq.48)) then ! 64x48

          print*,'water ice caps distribution for 64x48 resolution'
          longwatercaptag(1:65)   = .true. ! central cap - core
          longwatercaptag(75:85)  = .true. ! central cap 
          longwatercaptag(93:114) = .true. ! central cap
!---------------------   OUTLIERS  ----------------------------
          if (.true.) then
          longwatercaptag(136)    = .true. ! outlier, lat = 78.75
          longwatercaptag(138)    = .true. ! outlier, lat = 78.75
          longwatercaptag(140)    = .true. ! outlier, lat = 78.75
          longwatercaptag(142)    = .true. ! outlier, lat = 78.75
          longwatercaptag(161)    = .true. ! outlier, lat = 78.75
          longwatercaptag(163)    = .true. ! outlier, lat = 78.75
          longwatercaptag(165)    = .true. ! outlier, lat = 78.75
          longwatercaptag(183)    = .true. ! outlier, lat = 78.75
          longwatercaptag(185)    = .true. ! outlier, lat = 78.75
          longwatercaptag(187)    = .true. ! outlier, lat = 78.75
          longwatercaptag(189)    = .true. ! outlier, lat = 78.75
          longwatercaptag(191)    = .true. ! outlier, lat = 78.75
          longwatercaptag(193)    = .true. ! outlier, lat = 78.75
          longwatercaptag(194)    = .true. ! outlier, lat = 75
          longwatercaptag(203)    = .true. ! outlier, lat = 75
          longwatercaptag(207)    = .true. ! outlier, lat = 75
          longwatercaptag(244)    = .true. ! outlier, lat = 75
          longwatercaptag(246)    = .true. ! outlier, lat = 75
          longwatercaptag(250)    = .true. ! outlier, lat = 75
          longwatercaptag(252)    = .true. ! outlier, lat = 75
          longwatercaptag(254)    = .true. ! outlier, lat = 75
          longwatercaptag(256)    = .true. ! outlier, lat = 75
          endif
!--------------------------------------------------------------       


            
         else if (klon_glo .ne. 1) then
        
          print*,'No predefined ice location for this resolution :',
     &           nbp_lon,nbp_lat-1
          print*,'Please change icelocationmode in surfini.F'
          print*,'Or add some new definitions ...'
          call abort
          
         endif

         do ig=1,klon_glo
          if (longwatercaptag(ig)) watercaptag_glo(ig) = .true.
         enddo


        ELSE IF (icelocationmode .eq. 3) THEN
      
         print*,'Surfini: ice caps defined by lat and lon values'

         do ig=1,klon_glo
         
c-------- Towards olympia planitia water caps -----------
c-------------------------------------------------------- 

          if ( ( ( lati_glo(ig)*180./pi .ge. 77.  ) .and. ! cap #2
     .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
     .           ( long_glo(ig)*180./pi .ge. 110. ) .and.
     .           ( long_glo(ig)*180./pi .le. 181. ) )
     .         .or.

     .         ( ( lati_glo(ig)*180./pi .ge. 75.  ) .and. ! cap #4 (Korolev crater)
     .           ( lati_glo(ig)*180./pi .le. 76.  ) .and.
     .           ( long_glo(ig)*180./pi .ge. 150. ) .and.
     .           ( long_glo(ig)*180./pi .le. 168. ) )
     .         .or.
     .         ( ( lati_glo(ig)*180./pi .ge. 77 ) .and. ! cap #5
     .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
     .           ( long_glo(ig)*180./pi .ge. -150.) .and.
     .           ( long_glo(ig)*180./pi .le. -110.) ) )
     .         then
             
               if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
              !    watercaptag(ig)=.true.
                  alternate = 1
               else
                  alternate = 0
               endif !end if alternate = 0
               
          endif

c----------- Opposite olympia planitia water cap --------
c-------------------------------------------------------- 

          if ( ( ( lati_glo(ig)*180./pi     .ge.  80 ) .and.
     .         ( lati_glo(ig)*180./pi     .le.  84 ) )
     .         .and.
     .       ( ( long_glo(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
     .         ( long_glo(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
!     .     ( ( ( long_glo(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
!     .         ( long_glo(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
!     .       ( ( long_glo(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
!     .         ( long_glo(ig)*180./pi .le. -70. ) ) ) ) then  !!! 64x48
        !   watercaptag_glo(ig)=.true.
          endif


c -------------------- Central cap ----------------------
c-------------------------------------------------------- 

          if (abs(lati_glo(ig)*180./pi).gt.80)
     .          watercaptag_glo(ig)=.true.
           
c--------------------------------------------------------
c--------------------------------------------------------
         end do ! of (klon_glo)


        ELSE
      
         print*, 'In surfini.F, icelocationmode is ', icelocationmode
         print*, 'It should be 1, 2 or 3.'
         call abort 

        ENDIF ! of if (icelocation)
       
       
        ! print caps locations - useful for plots too 
        print*,'surfini: latitude | longitude | ig'
        do ig=1,klon_glo
          dryness_glo(ig) = icedryness

          if (watercaptag_glo(ig)) then
             print*,'surfini: ice water cap', lati_glo(ig)*180./pi,
     &              long_glo(ig)*180./pi, ig
          endif
        enddo
       
       endif !of if (is_master)
       
       ! Now scatter fields watercaptag and dryness from master to all
       ! (is just a plain copy in serial mode)
       call scatter(dryness_glo,dryness)
       call scatter(watercaptag_glo,watercaptag)
       
! end of #else of #ifdef MESOSCALE
       ENDIF ! (caps & water)
#endif       

c ===============================================================
c      INITIAL ALBEDO
c ===============================================================

         write(*,*)"surfini: water frost thickness",
     s     frost_albedo_threshold
         write(*,*)"surfini: water ice albedo:", albedo_h2o_ice
         write(*,*)"surfini: water ice TI:", inert_h2o_ice

c        To start with : Initial albedo = observed dataset 
c        -------------------------------------------------
         DO ig=1,ngrid
              psolaralb(ig,1)=albedodat(ig)
              psolaralb(ig,2)=albedodat(ig)
         END DO
         PRINT*,'surfini: minimum albedo without water caps',
     &          minval(albedodat)
         PRINT*,'surfini: maximum albedo without water caps',
     &          maxval(albedodat)

c        initial albedo if permanent H2O ice is present 
c        ------------------------------------------------
         IF ((water) .and. (caps)) THEN
           DO ig=1,ngrid
            IF (watercaptag(ig)) THEN
              psolaralb(ig,1) = albedo_h2o_ice
              psolaralb(ig,2) = albedo_h2o_ice
            ENDIF
           END DO
           PRINT*,'surfini: minimum albedo with water caps',
     &            minval(albedodat)
           PRINT*,'surfini: maximum albedo with water caps',
     &            maxval(albedodat)

         ENDIF

c      changing initial albedo if CO2 ice is present 
c      -------------------------------------------

       DO ig=1,ngrid
         IF (piceco2(ig) .GT. 0.) THEN
             IF(latitude(ig).LT. 0.) THEN
                icap=2 ! Southern hemisphere
             ELSE
                icap=1 ! Northern hemisphere
             ENDIF
             psolaralb(ig,1) = albedice(icap)
             psolaralb(ig,2) = albedice(icap)    
         END IF
       END DO

c      changing initial albedo if water ice frost is present
c      -------------------------------------------
       IF (water) THEN
          do iq=1,nqmx
c          if there is frost and surface albedo is set to albedo_h2o_ice
           if(noms(iq).eq."h2o_ice") then
             do ig=1,ngrid
             
              if ((watercaptag(ig).eqv..false.) 
     &     .and. (qsurf(ig,iq).lt.-frost_albedo_threshold)) then
              print*, ''
              print*, '!!! PROBLEM in SURFINI !!!!'
              print*, 'FOUND NEGATIVE SURFACE ICE VALUE WHERE
     & WATERCAPTAG IS FALSE'
              print*, ''
              print*, 'ig,qsurf,threshold' , 
     &         ig, qsurf(ig,iq), -frost_albedo_threshold
              print*, ''
              print*, '1) Check h2o_ice in startfi and ice
     & distribution in surfini'
              print*, '2) Use ini_h2osurf option in newstart'
              print*, ''
#ifndef MESOINI
              CALL ABORT
#else
              qsurf(ig,iq) = 0.
#endif
            endif
             
              if ((piceco2(ig) .eq. 0.).and.
     &          (qsurf(ig,iq).gt.frost_albedo_threshold)) then
                     psolaralb(ig,1) = albedo_h2o_ice
                     psolaralb(ig,2) = albedo_h2o_ice
c                     PRINT*,'surfini.F frost',
c     &                  latitude(ig)*180./pi, longitude(ig)*180./pi
               endif 
              enddo
           endif 
          end do
          PRINT*,'surfini: minimum albedo with frost and co2',
     &            minval(albedodat)
          PRINT*,'surfini: maximum albedo with frost and co2',
     &            maxval(albedodat)
       ELSE
         watercaptag(:) = .false.
       END IF ! OF IF(water)
         
      RETURN
      END