Ignore:
Timestamp:
Nov 29, 2023, 9:38:04 AM (19 months ago)
Author:
jbclement
Message:

Mars PCM:
Small "cosmetic" cleanings in some subroutines for the readability and "surfini.F" is transformed into "surfini_mod.F90" (explicit module interface + Fortran 90 format).
JBC

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
5 edited
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F

    r3130 r3142  
    9090
    9191      REAL,INTENT(INOUT) :: piceco2(ngrid,nslope) ! CO2 ice on the surface (kg.m-2)
    92       REAL,INTENT(INOUT) :: perennial_co2ice(ngrid,nslope) ! perennial CO2 ice on the surface (kg.m-2)
     92      REAL,INTENT(INOUT) :: perennial_co2ice(ngrid,nslope) ! Perennial CO2 ice on the surface (kg.m-2)
    9393      REAL,INTENT(INOUT) :: psolaralb(ngrid,2,nslope) ! albedo of the surface
    9494      REAL,INTENT(INOUT) :: pemisurf(ngrid,nslope) ! emissivity of the surface

  • trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90

    r3130 r3142  
    7979real, dimension(1)                  :: wstar = 0.    ! Thermals vertical velocity
    8080
    81 ! Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
     81! Physical and dynamical tendencies (e.g. m.s-2, K/s, Pa/s)
    8282real, dimension(nlayer)             :: du, dv, dtemp, dudyn, dvdyn, dtempdyn
    8383real, dimension(1)                  :: dpsurf

  • trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90

    r3132 r3142  
    738738endif !startphy_file
    739739
    740 
    741740if (paleoclimate) then
    742741  do iq=1,nq
    … …  
    744743   if (txt.eq."co2") igcm_co2_tmp = iq
    745744  enddo
     745
    746746  if (startphy_file) then
    747747! Depth of H2O lag
    … …  
    768768     lag_co2_ice(:,:) = -1.
    769769   endif
    770   else ! no startfiphyle
    771      h2o_ice_depth(:,:) = -1.
    772      lag_co2_ice(:,:) = -1.
    773      d_coef(:,:)= 4.e-4
    774   endif !startphy_file
    775 else
    776   h2o_ice_depth(:,:) = -1.
    777   lag_co2_ice(:,:) = -1.
    778   d_coef(:,:)= 4.e-4
    779 endif !paleoclimate
    780770
    781771! Perennial CO2 ice
    782 perennial_co2ice(:,:) = 0.
    783 if (paleoclimate) then
    784   if (startphy_file) then
     772    perennial_co2ice(:,:) = 0.
    785773    call get_field("perennial_co2ice",perennial_co2ice,found,indextime)
    786774    if (.not. found) then
    … …  
    794782        endif
    795783    endif ! not found
    796   else
     784  else ! no startfiphyle
     785    h2o_ice_depth(:,:) = -1.
     786    lag_co2_ice(:,:) = -1.
     787    d_coef(:,:)= 4.e-4
     788    perennial_co2ice(:,:) = 0.
    797789    if (abs(latitude(ngrid) - (-pi/2.)) < 1.e-5) then
    798790        do islope = 1,nslope
    … …  
    802794    endif
    803795  endif !startphy_file
    804 endif ! of if (paleoclimate)
     796else
     797   h2o_ice_depth(:,:) = -1.
     798   lag_co2_ice(:,:) = -1.
     799   d_coef(:,:)= 4.e-4
     800   perennial_co2ice(:,:) = 0.
     801endif !paleoclimate
    805802
    806803! close file:
    807 !
    808804if (startphy_file) call close_startphy
    809805

  • trunk/LMDZ.MARS/libf/phymars/phyredem.F90

    r3132 r3142  
    33implicit none
    44
     5!=======================================================================
    56contains
     7!=======================================================================
    68
    79subroutine physdem0(filename,lonfi,latfi,nsoil,ngrid,nlay,nq, &
    … …  
    911                         alb,ith,def_slope, &
    1012                         subslope_dist)
     13
    1114! create physics restart file and write time-independent variables
    1215  use comsoil_h, only: inertiedat, volcapa, mlayer
    … …  
    2528  use time_phylmdz_mod, only: daysec
    2629  use comslope_mod, ONLY: nslope
    27   USE paleoclimate_mod, ONLY: paleoclimate, h2o_ice_depth, lag_co2_ice 
     30  USE paleoclimate_mod, ONLY: paleoclimate, h2o_ice_depth, lag_co2_ice
     31
    2832  implicit none
    2933 
    … …  
    139143  ! Surface roughness length
    140144  call put_field("z0","Surface roughness length",z0)
    141  
    142   do ig=1,ngrid
    143     if(watercaptag(ig)) then
    144        watercaptag_tmp(ig)=1
    145     else
    146        watercaptag_tmp(ig)=0
    147     endif
    148   enddo
    149 
     145
     146  ! Water cap
     147  watercaptag_tmp = 0
     148  where (watercaptag) watercaptag_tmp = 1
    150149  call put_field("watercaptag","Infinite water reservoir",watercaptag_tmp)
    151150
    … …  
    155154
    156155  ! Paleoclimate outputs
    157   if(paleoclimate) then
     156  if (paleoclimate) then
    158157    call put_field("h2o_ice_depth","Depth to the fisrt H2O ice",h2o_ice_depth)
    159158    call put_field("lag_co2_ice","Depth of the CO2 lags",lag_co2_ice)
    160159  endif
     160
    161161  ! Close file
    162162  call close_restartphy
    163163 
    164164end subroutine physdem0
     165
     166!=======================================================================
    165167
    166168subroutine physdem1(filename,nsoil,ngrid,nlay,nq,nqsoil, &
    … …  
    227229 
    228230  ! Perennial CO2 ice layer
    229   if(paleoclimate) then
    230     call put_field("perennial_co2ice","CO2 ice cover",perennial_co2ice,time)
    231   endif
     231  if (paleoclimate) call put_field("perennial_co2ice","CO2 ice cover",perennial_co2ice,time)
    232232
    233233  ! Surface temperature

  • trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

    r3130 r3142  
    9191      use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
    9292      use check_fields_mod, only: check_physics_fields
     93      use surfini_mod, only: surfini
    9394#ifdef MESOSCALE
    9495      use comsoil_h, only: mlayer,layer

  • trunk/LMDZ.MARS/libf/phymars/surfini_mod.F90

    r3141 r3142  
     1MODULE surfini_mod
     2
     3implicit none
     4
     5!=======================================================================
     6contains
     7!=======================================================================
     8
    19      SUBROUTINE surfini(ngrid,nslope,qsurf)
    210
    … …  
    412      use netcdf
    513      use tracer_mod, only: nqmx, noms
    6       use geometry_mod, only: longitude, latitude, ! in radians
    7      &                     cell_area ! for watercaptag diagnosis
    8       use surfdat_h, only: watercaptag, frost_albedo_threshold,
    9      &                     albedo_h2o_cap, inert_h2o_ice, albedodat,
    10      &                     albedice, dryness
     14      use geometry_mod, only: longitude, latitude, & ! in radians
     15                          cell_area ! for watercaptag diagnosis
     16      use surfdat_h, only: watercaptag, frost_albedo_threshold, &
     17                          albedo_h2o_cap, inert_h2o_ice, albedodat, &
     18                          albedice, dryness
    1119#ifndef MESOSCALE
    1220      use mod_grid_phy_lmdz, only : klon_glo ! # of physics point on full grid
    … …  
    1624      use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
    1725      use datafile_mod, only: datadir
     26
    1827      IMPLICIT NONE
    19 c=======================================================================
    20 c
    21 c   creation des calottes pour l'etat initial
    22 c
    23 c=======================================================================
    24 c-----------------------------------------------------------------------
    25 c   Declarations:
    26 c   -------------
     28!=======================================================================
     29!
     30!   creation of caps for initial state
     31!
     32!=======================================================================
     33!-----------------------------------------------------------------------
     34!   Declarations:
     35!   -------------
    2736      include "callkeys.h"
    2837
    … …  
    3342      INTEGER ig,icap,iq,alternate,islope
    3443      REAL icedryness ! ice dryness
    35      
     44
    3645      ! longwatercaptag is watercaptag. Trick for some compilers
    3746      LOGICAL, DIMENSION(100000) :: longwatercaptag
    38      
     47
    3948! There are 4 different modes for ice distribution:
    4049! icelocationmode = 1 ---> based on data from surface.nc
    … …  
    4251! icelocationmode = 3 ---> based on logical relations for latitude and longitude
    4352! icelocationmode = 4 ---> predefined 64x48 but usable with every
    44 ! resolution, and easily adaptable for dynamico 
     53! resolution, and easily adaptable for dynamico
    4554! For visualisation : > /u/tnalmd/bin/watercaps gcm_txt_output_file
    4655      INTEGER,SAVE :: icelocationmode = 5
    47      
     56
    4857!$OMP THREADPRIVATE(icelocationmode)
    49        
    50        
     58
     59
    5160      !in case icelocationmode == 1
    5261      INTEGER i,j
    … …  
    5463      PARAMETER   (imd=360,jmd=180)
    5564      REAL        zdata(imd,jmd)
    56       REAL        zelat,zelon 
     65      REAL        zelat,zelon
    5766
    5867#ifndef MESOSCALE
    … …  
    6473
    6574      INTEGER   ierr,nid,nvarid
    66      
     75
    6776      REAL,SAVE :: min_icevalue = 500.
    68      
     77
    6978!$OMP THREADPRIVATE(min_icevalue)
    70      
     79
    7180      character(len=50) :: string = 'thermal'
    72      
     81
    7382      character (len=100) :: zedatafile
    7483
    … …  
    8291              qsurf(ig,iq,islope)=0.  !! on jette les inputs GCM
    8392                          !! on regle juste watercaptag
    84                           !! il faudrait garder les inputs GCM 
     93                          !! il faudrait garder les inputs GCM
    8594                          !! si elles sont consequentes
    8695            enddo
    8796         enddo
    88          if ( ( latitude(ig)*180./pi .gt. 70. ) .and.
    89      .        ( albedodat(ig) .ge. 0.26   ) )  then
     97         if ( ( latitude(ig)*180./pi .gt. 70. ) .and. ( albedodat(ig) .ge. 0.26   ) )  then
    9098                 write(*,*)"outlier ",ig
    9199                 watercaptag(ig)  = .true.
    92100                 dryness(ig)      = 1.
    93                  albedodat(ig)    = albedo_h2o_cap  !! pour output 
     101                 albedodat(ig)    = albedo_h2o_cap  !! pour output
    94102         else
    95103                 watercaptag(ig)  = .false.
    … …  
    118126#ifndef MESOSCALE
    119127
    120 c
    121 c=======================================================================
     128!
     129!=======================================================================
    122130! Initialize watercaptag (default is false)
    123131!      watercaptag_glo(:)=.false. !Already done in phyetat0 if needed
    124132
    125 c     water ice outliers
    126 c     ------------------------------------------
    127 
    128       IF (water) THEN 
    129      
    130 c Perennial H20 north cap defined by watercaptag=true (allows surface to be
    131 c hollowed by sublimation in vdifc).
    132 
    133 c We might not want albedodat to be modified because it is used to write
    134 c restart files. Instead, albedo is directly modified when needed (i.e.
    135 c if we have watercaptag and no co2 ice), below and in albedocaps.F90
    136 
    137 c       "Dryness coefficient" controlling the evaporation and
    138 c        sublimation from the ground water ice (close to 1)
    139 c        HERE, the goal is to correct for the fact
    140 c        that the simulated permanent water ice polar caps
    141 c        is larger than the actual cap and the atmospheric
    142 c        opacity not always realistic.
     133!     water ice outliers
     134!     ------------------------------------------
     135
     136      IF (water) THEN
     137
     138! Perennial H20 north cap defined by watercaptag=true (allows surface to be
     139! hollowed by sublimation in vdifc).
     140
     141! We might not want albedodat to be modified because it is used to write
     142! restart files. Instead, albedo is directly modified when needed (i.e.
     143! if we have watercaptag and no co2 ice), below and in albedocaps.F90
     144
     145!       "Dryness coefficient" controlling the evaporation and
     146!        sublimation from the ground water ice (close to 1)
     147!        HERE, the goal is to correct for the fact
     148!        that the simulated permanent water ice polar caps
     149!        is larger than the actual cap and the atmospheric
     150!        opacity not always realistic.
    143151
    144152         alternate = 0
    145          
     153
    146154         if (ngrid .ne. 1) then
    147155           IF (icelocationmode .ne. 5) then
    … …  
    150158           longwatercaptag(:) = .false.
    151159         endif
    152          
     160
    153161         write(*,*) "surfini: Ice dryness ?"
    154162         icedryness=1. ! default value
    … …  
    156164         write(*,*) "surfini: icedryness = ",icedryness
    157165         dryness (:) = icedryness
    158          
     166
    159167      ! To be able to run in parallel, we work on the full grid
    160168      ! and dispatch results afterwards
    … …  
    168176      if (is_master) then
    169177
    170         IF (ngrid .eq. 1) THEN ! special case for 1d --> do nothing
    171      
    172          print*, 'ngrid = 1, do no put ice caps in surfini.F'
    173 
    174         ELSE IF (icelocationmode .eq. 1) THEN
    175      
    176          print*,'Surfini: ice caps defined from surface.nc'
    177            
     178        IF (ngrid == 1) THEN ! special case for 1d --> do nothing
     179
     180         write(*,*) 'ngrid = 1, do no put ice caps in surfini.F'
     181
     182        ELSE
     183
     184        select case (icelocationmode)
     185            case(1) ! icelocationmode == 1
     186
     187         write(*,*)'Surfini: ice caps defined from surface.nc'
     188
    178189! This method detects ice as gridded value above min_icevalue in the field "string" from surface.nc
    179190! Typically, it is for thermal inertia above 500 tiu.
    … …  
    182193! (the approximation is that all points within the GCM latitude resolution have the same area).
    183194! 2. caps are placed to fill the GCM points with the most detected ice first.
    184      
    185 
    186            
     195
    187196         zedatafile = trim(datadir)
    188  
    189        
    190          ierr=nf90_open(trim(zedatafile)//'/surface.nc',
    191      &   NF90_NOWRITE,nid)
    192      
     197
     198         ierr=nf90_open(trim(zedatafile)//'/surface.nc',NF90_NOWRITE,nid)
     199
    193200         IF (ierr.NE.nf90_noerr) THEN
    194201       write(*,*)'Error : cannot open file surface.nc '
    … …  
    202209       call abort_physic("surfini","missing surface.nc file",1)
    203210         ENDIF
    204      
    205      
     211
     212
    206213         ierr=nf90_inq_varid(nid, string, nvarid)
    207214         if (ierr.ne.nf90_noerr) then
    … …  
    219226          call abort_physic("surfini","failed loading "//trim(string),1)
    220227         endif
    221  
    222                      
     228
     229
    223230         ierr=nf90_close(nid)
    224  
     231
    225232
    226233         nb_ice(:,1) = 1 ! default: there is no ice
    … …  
    230237         latice(:,2) = 0
    231238         lonice(:,2) = 0
    232          !print*,'jjm,iim',jjm,iim ! jjm =  nb lati , iim = nb longi
     239         !write(*,*)'jjm,iim',jjm,iim ! jjm =  nb lati , iim = nb longi
    233240
    234241         ! loop over the GCM grid - except for poles (ig=1 and ngrid)
    235242         do ig=2,klon_glo-1
    236      
    237           ! loop over the surface file grid     
     243
     244          ! loop over the surface file grid
    238245          do i=1,imd
    239246           do j=1,jmd
    240247             zelon = i - 180.
    241              zelat = 90. - j
    242             if ((abs(lati_glo(ig)*180./pi-zelat).le.
    243      &           90./real(nbp_lat-1)) .and.
    244      &          (abs(long_glo(ig)*180./pi-zelon).le.
    245      &           180./real(nbp_lon))) then
     248             zelat = 90. - j
     249            if ((abs(lati_glo(ig)*180./pi-zelat) .le. 90./real(nbp_lat-1)) .and. &
     250                (abs(long_glo(ig)*180./pi-zelon) .le. 180./real(nbp_lon))) then
    246251              ! count all points in that GCM grid point
    247252              nb_ice(ig,1) = nb_ice(ig,1) + 1
    248               if (zdata(i,j) > min_icevalue)
    249                  ! count all detected points in that GCM grid point
    250      &           nb_ice(ig,2) = nb_ice(ig,2) + 1
     253              ! count all detected points in that GCM grid point
     254              if (zdata(i,j) > min_icevalue) nb_ice(ig,2) = nb_ice(ig,2) + 1
    251255             endif
    252256           enddo
    253           enddo 
     257          enddo
    254258
    255259        ! projection of nb_ice on GCM lat and lon axes
    256           latice(1+(ig-2)/nbp_lon,:) =
    257      &     latice(1+(ig-2)/nbp_lon,:) + nb_ice(ig,:)
    258           lonice(1+mod(ig-2,nbp_lon),:) =
    259      &     lonice(1+mod(ig-2,nbp_lon),:) + nb_ice(ig,:) ! lonice is USELESS ...
     260          latice(1+(ig-2)/nbp_lon,:) = latice(1+(ig-2)/nbp_lon,:) + nb_ice(ig,:)
     261          lonice(1+mod(ig-2,nbp_lon),:) = lonice(1+mod(ig-2,nbp_lon),:) + nb_ice(ig,:) ! lonice is USELESS ...
    260262
    261263         enddo ! of do ig=2,klon_glo-1
    262      
    263 
    264      
     264
     265
     266
    265267         ! special case for poles
    266268         nb_ice(1,2)   = 1  ! ice prescribed on north pole
    … …  
    269271         latice(nbp_lat-1,:) = nb_ice(ngrid,:)
    270272         lonice(nbp_lon,:) = nb_ice(ngrid,:)
    271      
    272      
    273 !      print*, 'latice TOT', latice(:,1)
    274 !      print*, 'latice FOUND', latice(:,2)
    275 !      print*, 'lonice TOT', lonice(:,1)
    276 !      print*, 'lonice FOUND', lonice(:,2)
    277      
    278 !      print*, 'lat ratio', int(real(latice(:,2))/real(latice(:,1))*iim)
    279 !      print*, 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))*jjm)
    280      
    281 !      print*,''
    282 !      print*,'sum lat', sum(latice(:,1)), sum(lonice(:,1))
    283 !      print*,'sum lon', sum(latice(:,2)), sum(lonice(:,2))
    284      
    285    
     273
     274
     275!      write(*,*) 'latice TOT', latice(:,1)
     276!      write(*,*) 'latice FOUND', latice(:,2)
     277!      write(*,*) 'lonice TOT', lonice(:,1)
     278!      write(*,*) 'lonice FOUND', lonice(:,2)
     279
     280!      write(*,*) 'lat ratio', int(real(latice(:,2))/real(latice(:,1))*iim)
     281!      write(*,*) 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))*jjm)
     282
     283!      write(*,*)''
     284!      write(*,*)'sum lat', sum(latice(:,1)), sum(lonice(:,1))
     285!      write(*,*)'sum lon', sum(latice(:,2)), sum(lonice(:,2))
     286
     287
    286288         ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
    287289         do i=1,(nbp_lat-1)/2
    … …  
    290292          ! ratiolat is the ratio of area covered by ice within this GCM latitude range
    291293          ratiolat  = real(latice(i,2))/real(latice(i,1))
    292           !print*,'i',i,(i-1)*iim+2,i*iim+1
    293      
     294          !write(*,*)'i',i,(i-1)*iim+2,i*iim+1
     295
    294296          ! put ice caps while there is not enough ice,
    295297          ! as long as the threshold is above 20%
    … …  
    298300           ! loop over GCM longitudes
    299301           do j=1,nbp_lon
    300             ! if the detected ice ratio in the GCM grid point 
     302            ! if the detected ice ratio in the GCM grid point
    301303            ! is more than 'step', then add ice
    302             if (real(nb_ice((i-1)*nbp_lon+1+j,2))
    303      &        / real(nb_ice((i-1)*nbp_lon+1+j,1)) .ge. step) then
     304            if (real(nb_ice((i-1)*nbp_lon+1+j,2)) / real(nb_ice((i-1)*nbp_lon+1+j,1)) .ge. step) then
    304305                  watercaptag_glo((i-1)*nbp_lon+1+j) = .true.
    305306                  count = count + 1
    306307            endif
    307308           enddo ! of do j=1,nbp_lon
    308            !print*, 'step',step,count,ratiolat*nbp_lon
     309           !write(*,*) 'step',step,count,ratiolat*nbp_lon
    309310           step = step - 0.01
    310311          enddo ! of do while
    311           !print*, 'step',step,count,ratiolat*nbp_lon
     312          !write(*,*) 'step',step,count,ratiolat*nbp_lon
    312313
    313314         enddo ! of do i=1,jjm/2
    314            
    315 
    316         ELSE IF (icelocationmode .eq. 2) THEN
    317      
    318          print*,'Surfini: predefined ice caps'
    319      
     315
     316
     317            case(2) ! icelocationmode == 2
     318
     319         write(*,*)'Surfini: predefined ice caps'
     320
    320321         if ((nbp_lon.eq.32).and.((nbp_lat-1).eq.24)) then ! 32x24
    321            
    322           print*,'water ice caps distribution for 32x24 resolution'
     322
     323          write(*,*)'water ice caps distribution for 32x24 resolution'
    323324          longwatercaptag(1:9)    = .true. ! central cap - core
    324325          longwatercaptag(26:33)  = .true. ! central cap
    … …  
    333334         else if ((nbp_lon.eq.64).and.((nbp_lat-1).eq.48)) then ! 64x48
    334335
    335           print*,'water ice caps distribution for 64x48 resolution'
     336          write(*,*)'water ice caps distribution for 64x48 resolution'
    336337          longwatercaptag(1:65)   = .true. ! central cap - core
    337           longwatercaptag(75:85)  = .true. ! central cap 
     338          longwatercaptag(75:85)  = .true. ! central cap
    338339          longwatercaptag(93:114) = .true. ! central cap
    339340!---------------------   OUTLIERS  ----------------------------
    … …  
    362363          longwatercaptag(256)    = .true. ! outlier, lat = 75
    363364          endif
    364 !--------------------------------------------------------------       
    365 
    366 
    367            
     365!--------------------------------------------------------------
     366
     367
     368
    368369         else if (klon_glo .ne. 1) then
    369        
    370           print*,'No predefined ice location for this resolution :',
    371      &           nbp_lon,nbp_lat-1
    372           print*,'Please change icelocationmode in surfini.F'
    373           print*,'Or add some new definitions ...'
    374           call abort_physic("surfini",
    375      &         "no pre-definitions for this resolution",1)
    376          
     370
     371          write(*,*)'No predefined ice location for this resolution :',nbp_lon,nbp_lat-1
     372          write(*,*)'Please change icelocationmode in surfini.F'
     373          write(*,*)'Or add some new definitions ...'
     374          call abort_physic("surfini","no pre-definitions for this resolution",1)
     375
    377376         endif
    378377
    … …  
    382381
    383382
    384         ELSE IF (icelocationmode .eq. 3) THEN
    385      
    386          print*,'Surfini: ice caps defined by lat and lon values'
     383            case(3) ! icelocationmode == 3
     384
     385         write(*,*)'Surfini: ice caps defined by lat and lon values'
    387386
    388387         do ig=1,klon_glo
    389          
    390 c-------- Towards olympia planitia water caps -----------
    391 c--------------------------------------------------------
    392 
    393           if ( ( ( lati_glo(ig)*180./pi .ge. 77.  ) .and. ! cap #2
    394      .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
    395      .           ( long_glo(ig)*180./pi .ge. 110. ) .and.
    396      .           ( long_glo(ig)*180./pi .le. 181. ) )
    397      .         .or.
    398 
    399      .         ( ( lati_glo(ig)*180./pi .ge. 75.  ) .and. ! cap #4 (Korolev crater)
    400      .           ( lati_glo(ig)*180./pi .le. 76.  ) .and.
    401      .           ( long_glo(ig)*180./pi .ge. 150. ) .and.
    402      .           ( long_glo(ig)*180./pi .le. 168. ) )
    403      .         .or.
    404      .         ( ( lati_glo(ig)*180./pi .ge. 77 ) .and. ! cap #5
    405      .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
    406      .           ( long_glo(ig)*180./pi .ge. -150.) .and.
    407      .           ( long_glo(ig)*180./pi .le. -110.) ) )
    408      .         then
    409              
     388
     389!-------- Towards olympia planitia water caps -----------
     390!--------------------------------------------------------
     391
     392          if ( ( ( lati_glo(ig)*180./pi .ge. 77.  ) .and.  & ! cap #2
     393                 ( lati_glo(ig)*180./pi .le. 80.  ) .and.  &
     394                 ( long_glo(ig)*180./pi .ge. 110. ) .and.  &
     395                 ( long_glo(ig)*180./pi .le. 181. ) ) .or. &
     396               ( ( lati_glo(ig)*180./pi .ge. 75.  ) .and.  & ! cap #4 (Korolev crater)
     397                 ( lati_glo(ig)*180./pi .le. 76.  ) .and.  &
     398                 ( long_glo(ig)*180./pi .ge. 150. ) .and.  &
     399                 ( long_glo(ig)*180./pi .le. 168. ) ) .or. &
     400               ( ( lati_glo(ig)*180./pi .ge. 77 ) .and.    & ! cap #5
     401                 ( lati_glo(ig)*180./pi .le. 80.  ) .and.  &
     402                 ( long_glo(ig)*180./pi .ge. -150.) .and.  &
     403                 ( long_glo(ig)*180./pi .le. -110.) ) ) then
     404
    410405               if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
    411406              !    watercaptag(ig)=.true.
    … …  
    414409                  alternate = 0
    415410               endif !end if alternate = 0
    416                
     411
    417412          endif
    418413
    419 c----------- Opposite olympia planitia water cap --------
    420 c--------------------------------------------------------
    421 
    422           if ( ( ( lati_glo(ig)*180./pi     .ge.  80 ) .and.
    423      .         ( lati_glo(ig)*180./pi     .le.  84 ) )
    424      .         .and.
    425      .       ( ( long_glo(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
    426      .         ( long_glo(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
    427 !     .     ( ( ( long_glo(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
    428 !     .         ( long_glo(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
    429 !     .       ( ( long_glo(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
    430 !     .         ( long_glo(ig)*180./pi .le. -70. ) ) ) ) then  !!! 64x48
     414!----------- Opposite olympia planitia water cap --------
     415!--------------------------------------------------------
     416
     417          if ( ( ( lati_glo(ig)*180./pi .ge. 80 ) .and.   &
     418               ( lati_glo(ig)*180./pi   .le. 84 ) ) .and. &
     419             ( ( long_glo(ig)*180./pi .lt. -95. ) .or.    &  !!! 32x24
     420               ( long_glo(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
     421!           ( ( ( long_glo(ig)*180./pi .ge. -29. ) .and.  &   !!! 64x48
     422!               ( long_glo(ig)*180./pi .le.  90. ) ) .or. &   !!! 64x48
     423!             ( ( long_glo(ig)*180./pi .ge. -77. ) .and.  &   !!! 64x48
     424!               ( long_glo(ig)*180./pi .le. -70. ) ) ) ) then !!! 64x48
    431425        !   watercaptag_glo(ig)=.true.
    432426          endif
    433427
    434428
    435 c -------------------- Central cap ----------------------
    436 c--------------------------------------------------------
    437 
    438           if (abs(lati_glo(ig)*180./pi).gt.80)
    439      .          watercaptag_glo(ig)=.true.
    440            
    441 c--------------------------------------------------------
    442 c--------------------------------------------------------
     429! -------------------- Central cap ----------------------
     430!--------------------------------------------------------
     431
     432          if (abs(lati_glo(ig)*180./pi).gt.80) watercaptag_glo(ig)=.true.
     433
     434!--------------------------------------------------------
     435!--------------------------------------------------------
    443436         end do ! of (klon_glo)
    444437
    445         ELSE IF (icelocationmode .eq. 4) THEN
    446      
    447          print*,'icelocationmode = 4'
    448          print*,'Surfini: ice caps defined using manual 64x48 settings'
    449          print*,'(although, it should work with any resolution)'
    450          call locate_watercaptag(klon_glo,lati_glo,
    451      &            long_glo,watercaptag_glo)
    452 
    453 !         print*,'watercaptag_glo(:), ',watercaptag_glo(:)
    454 
    455         ELSE IF (icelocationmode .eq. 5) THEN
    456      
    457          print*,'icelocationmode = 5'
    458          print*,'Surfini: ice caps defined using startfi.nc data'
     438            case(4) ! icelocationmode == 4
     439
     440         write(*,*)'icelocationmode = 4'
     441         write(*,*)'Surfini: ice caps defined using manual 64x48 settings'
     442         write(*,*)'(although, it should work with any resolution)'
     443         call locate_watercaptag(klon_glo,lati_glo,long_glo,watercaptag_glo)
     444
     445!         write(*,*)'watercaptag_glo(:), ',watercaptag_glo(:)
     446
     447            case(5) ! icelocationmode == 5
     448
     449         write(*,*)'icelocationmode = 5'
     450         write(*,*)'Surfini: ice caps defined using startfi.nc data'
    459451         do ig=1,ngrid
    460452           if(any(watercaptag_glo)) then
    461453           else
    462               call locate_watercaptag(klon_glo,lati_glo,
    463      &            long_glo,watercaptag_glo)
     454              call locate_watercaptag(klon_glo,lati_glo,long_glo,watercaptag_glo)
    464455           endif
    465456         enddo
    466457
    467 !         print*,'watercaptag_glo(:), ',watercaptag_glo(:)
    468 
    469         ELSE
    470      
    471          print*, 'In surfini.F, icelocationmode is ', icelocationmode
    472          print*, 'It should be 1, 2, 3, 4 or 5 (default is 5)'
     458!         write(*,*)'watercaptag_glo(:), ',watercaptag_glo(:)
     459
     460            case default
     461
     462         write(*,*) 'In surfini.F, icelocationmode is ', icelocationmode
     463         write(*,*) 'It should be 1, 2, 3, 4 or 5 (default is 5)'
    473464         call abort_physic("surfini","wrong icelocationmode",1)
    474465
    475         ENDIF ! of if (icelocation)
    476        
    477        
    478         ! print caps locations - useful for plots too
    479         print*,'surfini: latitude | longitude | ig'
     466        end select
     467
     468        ENDIF ! of if (ngrid)
     469
     470
     471        ! print caps locations - useful for plots too
     472        write(*,*)'surfini: latitude | longitude | ig'
    480473        do ig=1,klon_glo
    481474          dryness_glo(ig) = icedryness
    482475
    483476          if (watercaptag_glo(ig)) then
    484              print*,'surfini: ice water cap', lati_glo(ig)*180./pi,
    485      &              long_glo(ig)*180./pi, ig
    486 !             write(1,*),ig, lati_glo(ig)*180./pi,
    487 !     &              cell_area(ig)
    488 !             write(2,*), lati_glo(ig)*180./pi,
    489 !     &              long_glo(ig)*180./pi,cell_area(ig)
    490 !             write(3,*), ig, lati_glo(ig)*180./pi,
    491 !     &              long_glo(ig)*180./pi,cell_area(ig)
     477             write(*,*)'surfini: ice water cap', lati_glo(ig)*180./pi,long_glo(ig)*180./pi, ig
     478!             write(1,*),ig, lati_glo(ig)*180./pi,cell_area(ig)
     479!             write(2,*), lati_glo(ig)*180./pi,long_glo(ig)*180./pi,cell_area(ig)
     480!             write(3,*), ig, lati_glo(ig)*180./pi,long_glo(ig)*180./pi,cell_area(ig)
    492481          endif
    493482        enddo
    494        
     483
    495484       endif !of if (is_master)
    496        
     485
    497486       if (ngrid.gt.1) then
    498487        ! Now scatter fields watercaptag and dryness from master to all
    … …  
    505494       ENDIF ! water
    506495! end of #else of #ifndef MESOSCALE
    507 #endif       
    508 !      END SUBROUTINE surfini(ngrid,piceco2,qsurf)
    509       END !SUBROUTINE surfini(ngrid,piceco2,qsurf)
    510 
    511       SUBROUTINE locate_watercaptag(klon_glo,lati_glo,
    512      &            long_glo,watercaptag_glo)
     496#endif
     497      END SUBROUTINE surfini
     498
     499!=======================================================================
     500
     501      SUBROUTINE locate_watercaptag(klon_glo,lati_glo,long_glo,watercaptag_glo)
    513502
    514503      USE comcstfi_h, ONLY: pi
    … …  
    527516! coordinates in latitude and longitude are written below. With this
    528517! routine, we check if the grid cell center is in between any of those
    529 ! points. If so, watercaptag = true.
    530 
    531 
    532 
    533 
    534       latedge(:,1)=(/
    535      & 88.125, 84.375, 84.375, 84.375, 84.375, 84.375,84.375, 84.375,
    536      & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
    537      & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
    538      & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
    539      & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
    540      & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
    541      & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
    542      & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
    543      & 84.375, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
    544      & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
    545      & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
    546      & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
    547      & 80.625, 80.625, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875,
    548      & 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 73.125,
    549      & 73.125, 73.125, 73.125, 73.125, 73.125, 73.125, 73.125, 73.125/)
    550 
    551 
    552       latedge(:,2)=(/
    553      & 90. , 88.125, 88.125, 88.125, 88.125, 88.125,88.125, 88.125,   
    554      & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
    555      & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
    556      & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
    557      & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
    558      & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
    559      & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
    560      & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
    561      & 88.125, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
    562      & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
    563      & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
    564      & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
    565      & 84.375, 84.375, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
    566      & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 76.875,
    567      & 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875/)
    568 
    569 
    570       lonedge(:,1)=(/
    571      &-180.    , -180.    , -177.1875, -171.5625,-165.9375, -160.3125,
    572      &-154.6875, -149.0625, -143.4375, -137.8125, -132.1875,-126.5625,
    573      &-120.9375, -115.3125, -109.6875, -104.0625,  -98.4375, -92.8125,
    574      & -87.1875,  -81.5625,  -75.9375,  -70.3125,  -64.6875, -59.0625,
    575      & -53.4375,  -47.8125,  -42.1875,  -36.5625,  -30.9375, -25.3125,
    576      & -19.6875,  -14.0625,   -8.4375,   -2.8125,    2.8125,   8.4375,
    577      &  14.0625,   19.6875,   25.3125,   30.9375,   36.5625,  42.1875,
    578      &  47.8125,   53.4375,   59.0625,   64.6875,   70.3125,  75.9375,
    579      &  81.5625,   87.1875,   92.8125,   98.4375,  104.0625, 109.6875,
    580      & 115.3125,  120.9375,  126.5625,  132.1875,  137.8125, 143.4375,
    581      & 149.0625,  154.6875,  160.3125,  165.9375,  171.5625,-132.1875,
    582      &-126.5625, -120.9375, -115.3125, -109.6875, -104.0625, -98.4375,
    583      & -92.8125,  -87.1875,  -81.5625,  -75.9375,  -30.9375, -25.3125,
    584      & -19.6875,  -14.0625,   -8.4375,   -2.8125,    2.8125,   8.4375,
    585      &  14.0625,   19.6875,   25.3125,   30.9375,   36.5625,  42.1875,
    586      &  47.8125,   53.4375,   59.0625,   64.6875,   70.3125,  75.9375,
    587      &  81.5625,   87.1875, -149.0625, -137.8125, -126.5625,-115.3125,
    588      &  -8.4375,    2.8125,   14.0625,  115.3125,  126.5625, 137.8125,
    589      & 149.0625,  160.3125,  171.5625, -180.    , -132.1875,-109.6875,
    590      &  98.4375,  109.6875,  132.1875,  143.4375,  154.6875,165.9375/)
    591 
    592       lonedge(:,2)=(/
    593      & 180.    , -180.    , -171.5625, -165.9375,-160.3125, -154.6875,
    594      &-149.0625,-143.4375, -137.8125, -132.1875, -126.5625, -120.9375,
    595      &-115.3125,-109.6875, -104.0625,  -98.4375,  -92.8125,  -87.1875,
    596      & -81.5625, -75.9375,  -70.3125,  -64.6875,  -59.0625,  -53.4375,
    597      & -47.8125, -42.1875,  -36.5625,  -30.9375,  -25.3125,  -19.6875,
    598      & -14.0625,  -8.4375,   -2.8125,    2.8125,    8.4375,   14.0625,
    599      &  19.6875,  25.3125,   30.9375,   36.5625,   42.1875,   47.8125,
    600      &  53.4375,  59.0625,   64.6875,   70.3125,   75.9375,   81.5625,
    601      &  87.1875,  92.8125,   98.4375,  104.0625,  109.6875,  115.3125,
    602      & 120.9375, 126.5625,  132.1875,  137.8125,  143.4375,  149.0625,
    603      & 154.6875, 160.3125,  165.9375,  171.5625,  177.1875, -126.5625,
    604      &-120.9375,-115.3125, -109.6875, -104.0625,  -98.4375,  -92.8125,
    605      & -87.1875, -81.5625,  -75.9375,  -70.3125,  -25.3125,  -19.6875,
    606      & -14.0625,  -8.4375,   -2.8125,    2.8125,    8.4375,   14.0625,
    607      &  19.6875,  25.3125,   30.9375,   36.5625,   42.1875,   47.8125,
    608      &  53.4375,  59.0625,   64.6875,   70.3125,   75.9375,   81.5625,
    609      &  87.1875,  92.8125, -143.4375, -132.1875, -120.9375, -109.6875,
    610      &  -2.8125,   8.4375,   19.6875,  120.9375,  132.1875,  143.4375,
    611      & 154.6875, 165.9375,  177.1875, -177.1875, -126.5625, -104.0625,
    612      & 104.0625, 115.3125,  137.8125,  149.0625,  160.3125,171.5625/)
    613 
     518! points. If so, watercaptag = true.
     519
     520      latedge(:,1)=(/                                                 &
     521      88.125, 84.375, 84.375, 84.375, 84.375, 84.375,84.375, 84.375,  &
     522      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
     523      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
     524      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
     525      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
     526      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
     527      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
     528      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
     529      84.375, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, &
     530      80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, &
     531      80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, &
     532      80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, &
     533      80.625, 80.625, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, &
     534      76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 73.125, &
     535      73.125, 73.125, 73.125, 73.125, 73.125, 73.125, 73.125, 73.125/)
     536
     537      latedge(:,2)=(/                                                 &
     538      90. , 88.125, 88.125, 88.125, 88.125, 88.125,88.125, 88.125,    &
     539      88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, &
     540      88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, &
     541      88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, &
     542      88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, &
     543      88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, &
     544      88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, &
     545      88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, &
     546      88.125, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
     547      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
     548      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
     549      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
     550      84.375, 84.375, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, &
     551      80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 76.875, &
     552      76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875/)
     553
     554      lonedge(:,1)=(/                                                 &
     555     -180.    , -180.    , -177.1875, -171.5625,-165.9375, -160.3125, &
     556     -154.6875, -149.0625, -143.4375, -137.8125, -132.1875,-126.5625, &
     557     -120.9375, -115.3125, -109.6875, -104.0625,  -98.4375, -92.8125, &
     558      -87.1875,  -81.5625,  -75.9375,  -70.3125,  -64.6875, -59.0625, &
     559      -53.4375,  -47.8125,  -42.1875,  -36.5625,  -30.9375, -25.3125, &
     560      -19.6875,  -14.0625,   -8.4375,   -2.8125,    2.8125,   8.4375, &
     561       14.0625,   19.6875,   25.3125,   30.9375,   36.5625,  42.1875, &
     562       47.8125,   53.4375,   59.0625,   64.6875,   70.3125,  75.9375, &
     563       81.5625,   87.1875,   92.8125,   98.4375,  104.0625, 109.6875, &
     564      115.3125,  120.9375,  126.5625,  132.1875,  137.8125, 143.4375, &
     565      149.0625,  154.6875,  160.3125,  165.9375,  171.5625,-132.1875, &
     566     -126.5625, -120.9375, -115.3125, -109.6875, -104.0625, -98.4375, &
     567      -92.8125,  -87.1875,  -81.5625,  -75.9375,  -30.9375, -25.3125, &
     568      -19.6875,  -14.0625,   -8.4375,   -2.8125,    2.8125,   8.4375, &
     569       14.0625,   19.6875,   25.3125,   30.9375,   36.5625,  42.1875, &
     570       47.8125,   53.4375,   59.0625,   64.6875,   70.3125,  75.9375, &
     571       81.5625,   87.1875, -149.0625, -137.8125, -126.5625,-115.3125, &
     572       -8.4375,    2.8125,   14.0625,  115.3125,  126.5625, 137.8125, &
     573      149.0625,  160.3125,  171.5625, -180.    , -132.1875,-109.6875, &
     574       98.4375,  109.6875,  132.1875,  143.4375,  154.6875,165.9375/)
     575
     576      lonedge(:,2)=(/                                                  &
     577       180.    , -180.    , -171.5625, -165.9375,-160.3125, -154.6875, &
     578      -149.0625,-143.4375, -137.8125, -132.1875, -126.5625, -120.9375, &
     579      -115.3125,-109.6875, -104.0625,  -98.4375,  -92.8125,  -87.1875, &
     580       -81.5625, -75.9375,  -70.3125,  -64.6875,  -59.0625,  -53.4375, &
     581       -47.8125, -42.1875,  -36.5625,  -30.9375,  -25.3125,  -19.6875, &
     582       -14.0625,  -8.4375,   -2.8125,    2.8125,    8.4375,   14.0625, &
     583        19.6875,  25.3125,   30.9375,   36.5625,   42.1875,   47.8125, &
     584        53.4375,  59.0625,   64.6875,   70.3125,   75.9375,   81.5625, &
     585        87.1875,  92.8125,   98.4375,  104.0625,  109.6875,  115.3125, &
     586       120.9375, 126.5625,  132.1875,  137.8125,  143.4375,  149.0625, &
     587       154.6875, 160.3125,  165.9375,  171.5625,  177.1875, -126.5625, &
     588      -120.9375,-115.3125, -109.6875, -104.0625,  -98.4375,  -92.8125, &
     589       -87.1875, -81.5625,  -75.9375,  -70.3125,  -25.3125,  -19.6875, &
     590       -14.0625,  -8.4375,   -2.8125,    2.8125,    8.4375,   14.0625, &
     591        19.6875,  25.3125,   30.9375,   36.5625,   42.1875,   47.8125, &
     592        53.4375,  59.0625,   64.6875,   70.3125,   75.9375,   81.5625, &
     593        87.1875,  92.8125, -143.4375, -132.1875, -120.9375, -109.6875, &
     594        -2.8125,   8.4375,   19.6875,  120.9375,  132.1875,  143.4375, &
     595       154.6875, 165.9375,  177.1875, -177.1875, -126.5625, -104.0625, &
     596       104.0625, 115.3125,  137.8125,  149.0625,  160.3125,171.5625/)
    614597
    615598      watercaptag_glo(:) = .false.
    616599      DO ig=1, klon_glo
    617600        DO i=1, 120
    618            if ((long_glo(ig)*180./pi.ge.lonedge(i,1))
    619      &         .and.(long_glo(ig)*180./pi.le.lonedge(i,2))
    620      &         .and.(lati_glo(ig)*180./pi.ge.latedge(i,1))
    621      &         .and.(lati_glo(ig)*180./pi.le.latedge(i,2))) then
     601           if ((long_glo(ig)*180./pi.ge.lonedge(i,1))      &
     602               .and.(long_glo(ig)*180./pi.le.lonedge(i,2)) &
     603               .and.(lati_glo(ig)*180./pi.ge.latedge(i,1)) &
     604               .and.(lati_glo(ig)*180./pi.le.latedge(i,2))) then
    622605             watercaptag_glo(ig) = .true.
    623606           endif
    … …  
    626609
    627610      END SUBROUTINE locate_watercaptag
     611
     612END MODULE surfini_mod
Note: See TracChangeset for help on using the changeset viewer.