Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/create_limit.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/create_limit.F	(revision 136)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/create_limit.F	(revision 173)
@@ -1,3 +1,6 @@
+c $Header$
       PROGRAM create_limit
+      USE startvar
+      USE ioipsl
       IMPLICIT none
 c
@@ -24,19 +27,19 @@
 #include "comgeom2.h"
 #include "comconst.h"
-c
+#include "dimphy.h"
+#include "indicesol.h"
 c-----------------------------------------------------------------------
-      INTEGER KIDIA, KFDIA, KLON, KLEV
-      PARAMETER (KIDIA=1,KFDIA=iim*(jjm-1)+2,
-     .           KLON=KFDIA-KIDIA+1,KLEV=llm)
-c-----------------------------------------------------------------------
-      REAL phy_nat(klon,360), phy_nat0(klon)
+      REAL phy_nat(klon,360)
+      real phy_nat0(klon)
       REAL phy_alb(klon,360)
       REAL phy_sst(klon,360)
       REAL phy_bil(klon,360)
       REAL phy_rug(klon,360)
-      REAL phy_ice(klon,360)
+      REAL phy_ice(klon)
 CPB
-      REAL phy_icet(klon,360)
-      REAL phy_oce(klon,360)
+c      REAL phy_icet(klon,360)
+c      REAL phy_oce(klon,360)
+      real pctsrf_t(klon,nbsrf,360)
+      real pctsrf(klon,nbsrf)
       REAL verif
 c
@@ -89,7 +92,25 @@
       INTEGER id_FOCE, id_FSIC, id_FTER, id_FLIC
 
-      INTEGER i, j, k, l
+      INTEGER i, j, k, l, ji
+c declarations pour lecture glace de mer
+      INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
+      INTEGER :: itaul(1), fid
+      REAL :: lev(1), date, dt
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic
+      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_lic, dlat_lic
+      REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
+      REAL :: flic_tmp(iip1, jjp1)
+      REAL :: champint(iim, jjp1)
 c Diverses variables locales
       REAL time
+! pour la lecture du fichier masque ocean
+      integer :: nid_o2a
+      logical :: couple = .false.
+      INTEGER :: iml_omask, jml_omask
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask
+      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_omask, dlat_omask
+      REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp
+      real, dimension(klon) :: ocemask_fi
+
 
       INTEGER          longcles
@@ -99,4 +120,5 @@
       INTEGER ncid,varid,ndimid(4),dimid
       character*30 namedim
+      CHARACTER*80 :: varname
 
 c initialisations:
@@ -122,64 +144,18 @@
 C Traitement du relief au sol
 c
-      write(*,*) 'Traitement du relief au sol pour fabriquer masque'
-      ierr = NF_OPEN('Relief.nc', NF_NOWRITE, ncid)
-
-      if (ierr.ne.0) then
-        print *, NF_STRERROR(ierr)
-        STOP
-      ENDIF
-
-      ierr = NF_INQ_VARID(ncid,'RELIEF',varid)
-      if (ierr.ne.0) then
-        print *, NF_STRERROR(ierr)
-        STOP
-      ENDIF
-      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
-      if (ierr.ne.0) then
-        print *, NF_STRERROR(ierr)
-        STOP
-      ENDIF
-      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
-      if (ierr.ne.0) then
-        print *, NF_STRERROR(ierr)
-        STOP
-      ENDIF
-      print*,'variable ', namedim, 'dimension ', imdep
-      ierr = NF_INQ_VARID(ncid,namedim,dimid)
-      if (ierr.ne.0) then
-        print *, NF_STRERROR(ierr)
-        STOP
-      ENDIF
-      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_msk)
-      if (ierr.ne.0) then
-        print *, NF_STRERROR(ierr)
-        STOP
-      ENDIF
-      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
-      if (ierr.ne.0) then
-        print *, NF_STRERROR(ierr)
-        STOP
-      ENDIF
-      print*,'variable ', namedim, 'dimension ', jmdep
-      ierr = NF_INQ_VARID(ncid,namedim,dimid)
-      if (ierr.ne.0) then
-        print *, NF_STRERROR(ierr)
-        STOP
-      ENDIF
-      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_msk)
-      if (ierr.ne.0) then
-        print *, NF_STRERROR(ierr)
-        STOP
-      ENDIF
-      ierr = NF_GET_VAR_REAL(ncid,varid,champ_msk)
-      if (ierr.ne.0) then
-        print *, NF_STRERROR(ierr)
-        STOP
-      ENDIF
-
-c
-      CALL mask_c_o(imdep, jmdep, dlon_msk, dlat_msk,champ_msk,
-     .             iim, jjp1, rlonv, rlatu, champint)
-      CALL gr_int_dyn(champint, masque, iim, jjp1)
+      write(*,*) 'Fabrication masque'
+
+      varname = 'masque'
+      masque(:,:) = 0.0
+      CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0)
+      pctsrf=0.
+      varname = 'zmasq'
+      zmasq(:) = 0.
+      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmasq,0.0)
+      WHERE (zmasq(1 : klon) .LE. EPSFRA)
+          zmasq(1 : klon) = 0.
+      END WHERE 
+!      WRITE(*,*)zmasq
+
       IF ( fracterre ) THEN
           DO i = 1, iim
@@ -195,9 +171,126 @@
       DO i = 1, iim
       DO j = 1, jjp1
-         mask(i,j) = champint(i,j)
+         mask(i,j) = masque(i,j)
       ENDDO
       ENDDO
       CALL gr_dyn_fi(1, iip1, jjp1, klon, masque, phy_nat0)
-      ierr = NF_CLOSE(ncid)
+C
+C En cas de simulation couplee, lecture du masque ocean issu du modele ocean
+C utilise pour calculer les poids et pour assurer l'adequation entre les
+C fractions d'ocean vu par l'atmosphere et l'ocean
+C
+
+      write(*,*)'Essai de lecture masque ocean'
+      iret = nf_open("o2a.nc", NF_NOWRITE, nid_o2a)
+      if (iret .ne. 0) then
+        write(*,*)'ATTENTION!! pas de fichier o2a.nc trouve'
+        write(*,*)'Run force'
+      else
+        couple = .true.
+        iret = nf_close(nid_o2a)
+        call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp
+     $    , nid_o2a)
+        if (iml_omask /= iim .or. jml_omask /= jjp1) then
+          write(*,*)'Dimensions non compatibles pour masque ocean'
+          write(*,*)'iim = ',iim,' iml_omask = ',iml_omask
+          write(*,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
+          stop
+        endif
+        ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret)
+        ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret)
+        ALLOCATE(dlon_omask(iml_omask), stat=iret)
+        ALLOCATE(dlat_omask(jml_omask), stat=iret)
+        ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret)
+        ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret)
+        CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp
+     $    , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
+        CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp, 
+     $      ttm_tmp, 1, 1, ocetmp)
+        CALL flinclo(fid)
+        dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1)
+        dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask)
+        ocemask = ocetmp
+        if (dlat_omask(1) < dlat_omask(jml_omask)) then
+          do j = 1, jml_omask
+            ocemask(:,j) = ocetmp(:,jml_omask-j+1)
+          enddo
+        endif 
+C
+C passage masque ocean a la grille physique
+C
+        ocemask_fi(1) = ocemask(1,1)
+        do j = 2, jjm
+          do i = 1, iim
+            ocemask_fi((j-2)*iim + i + 1) = ocemask(i,j)
+          enddo
+        enddo
+        ocemask_fi(klon) = ocemask(1,jjp1)
+        zmasq = 1. - ocemask_fi
+      endif
+
+
+C
+C lecture du fichier glace de terre pour fixer la fraction de terre 
+C et de glace de terre
+C
+      CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp
+     $    , fid)
+      ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret)
+      ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret)
+      ALLOCATE(dlon_lic(iml_lic), stat=iret)
+      ALLOCATE(dlat_lic(jml_lic), stat=iret)
+      ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret)
+      CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp
+     $    , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
+      CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp
+     $    , 1, 1, fraclic)
+      CALL flinclo(fid)
+C
+C interpolation sur la grille T du modele
+C
+      WRITE(*,*) 'dimensions de landice iml_lic, jml_lic : ', 
+     $    iml_lic, jml_lic
+c
+C sil les coordonnees sont en degres, on les transforme
+C
+      IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) )  THEN
+          lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180.
+      ENDIF 
+      IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN 
+          lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180.
+      ENDIF 
+
+      dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1)
+      dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic) 
+C
+      CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic
+     $    ,iim, jjp1,
+     $    rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1))
+c$$$      flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1)
+      flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1)
+C
+C passage sur la grille physique
+C
+      CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp,
+     $    pctsrf(1:klon, is_lic))
+C adequation avec le maque terre/mer
+      WHERE (pctsrf(1 : klon, is_lic) .LE. EPSFRA ) 
+          pctsrf(1 : klon, is_lic) = 0. 
+      END WHERE
+      WHERE (zmasq( 1 : klon) .LE. EPSFRA) 
+          pctsrf(1 : klon, is_lic) = 0.
+      END WHERE 
+      pctsrf(1 : klon, is_ter) = zmasq(1 : klon)
+      DO ji = 1, klon
+        IF (zmasq(ji) .GT. EPSFRA) THEN 
+            IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN
+                pctsrf(ji, is_lic) = zmasq(ji)
+                pctsrf(ji, is_ter) = 0.
+            ELSE 
+                pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
+            ENDIF 
+        ENDIF 
+      END DO 
+
 c
 c
@@ -295,4 +388,5 @@
          ENDDO
       ENDDO
+c      write(70,*)champtime
 c
       DO l = 1, lmdep
@@ -453,112 +547,77 @@
       ENDDO
 c
-      WRITE(*,*) 'phy_nat'
+c      WRITE(*,*) 'phy_nat'
 c     WRITE(*,'(72f4.1)') phy_nat0(1:klon)
 c
       DO k = 1, 360
         CALL gr_dyn_fi(1, iip1, jjp1, klon,
-     .      champan(1,1,k), phy_ice(1,k))
+     .      champan(1,1,k), phy_ice)
         IF ( newlmt) THEN
 
 CPB  en attendant de mettre fraction de terre
 c
-            WHERE(phy_ice(1:klon, k) .GT. 1.) phy_ice(1 : klon, k) = 1.
-            WHERE(phy_ice(1:klon, k) .LT. 0.) phy_ice(1 : klon, k) = 0.
-            WRITE(*,*) 'phy_ice : ', k
-c           WRITE(*,'(72f4.1)') phy_ice(1 : klon, k)
+          WHERE(phy_ice(1:klon) .GT. 1.) phy_ice(1 : klon) = 1.
+          WHERE(phy_ice(1:klon) .LT. EPSFRA) phy_ice(1 : klon) = 0.
 c 
-            IF (fracterre ) THEN
-                WRITE(*,*) 'passe dans cas fracterre' 
-                DO i = 1, klon
-                  phy_nat(i,k) = phy_nat0(i)
-                  IF (phy_nat0(i) .GE. 0.5 ) THEN
-                      IF(phy_ice(i,k) .GE. 1.e-5) THEN 
-                          IF ( phy_ice(i,k) .LE. phy_nat(i,k)) THEN
-                              phy_oce(i,k) = 1. - phy_nat(i,k)
-                              phy_icet(i,k) = phy_ice(i,k)
-                              phy_ice(i,k) = 0.
-                              phy_nat(i,k)= phy_nat(i,k)- phy_icet(i,k)
-                          ELSE 
-                              phy_oce(i,k) = 1. - phy_ice(i,k)
-                              phy_icet(i,k) = phy_nat(i,k)
-                              phy_nat(i,k) = 0.
-                              phy_ice(i,k) = phy_ice(i,k) 
-     $                            - phy_icet(i,k)
-                          ENDIF 
-                      ELSE
-                          phy_icet(i,k) = 0.
-                          phy_ice(i,k) = 0.
-                          phy_oce(i,k) = 1. - phy_nat(i,k)
-                      ENDIF 
-                  ELSE
-                      phy_oce(i,k) = 1. - phy_nat(i,k)
-                      IF(phy_ice(i,k) .GE. 1.e-5) THEN 
-                          IF( phy_ice(i,k) .LE. phy_oce(i,k) ) THEN 
-                              phy_icet(i,k) = 0.
-                              phy_oce(i,k) = phy_oce(i,k) - phy_ice(i,k)
-                          ELSE
-                              phy_icet(i,k)=phy_ice(i,k) - phy_oce(i,k)
-                              phy_ice(i,k) = phy_oce(i,k)
-                              phy_oce(i,k) = 0.
-                              phy_nat(i,k) = phy_nat(i,k)-phy_icet(i,k)
-                          ENDIF 
-                      ELSE
-                          phy_icet(i,k) = 0.
-                          phy_ice(i,k) = 0.
-                      ENDIF 
-                  ENDIF
-                  verif = phy_nat(i,k) + phy_icet(i,k)+ phy_ice(i,k) 
-     $                + phy_oce(i,k)
-                  IF ( verif .LT. 1. - 1.e-5 .OR. 
-     $                verif .GT. 1 + 1.e-5) THEN  
-                      WRITE(*,*) 'pb sous maille au point : i,k,verif '
+          IF (fracterre ) THEN
+c            WRITE(*,*) 'passe dans cas fracterre' 
+            pctsrf_t(:,is_ter,k) = pctsrf(:,is_ter)
+            pctsrf_t(:,is_lic,k) = pctsrf(:,is_lic)
+            DO i = 1, klon
+              pctsrf_t(i,is_sic,k) = (1. - pctsrf_t(i,is_lic,k) -
+     .                               pctsrf_t(i,is_ter,k)) * phy_ice(i)
+              pctsrf_t(i,is_oce,k) = 1. - pctsrf_t(i,is_lic,k) -
+     .                      pctsrf_t(i,is_ter,k) - pctsrf_t(i,is_sic,k)  
+              if (pctsrf_t(i,is_oce,k) .lt. 0.) then
+                      WRITE(*,*) 'pb sous maille au point : i,k '
+     $                    , i,k,pctsrf_t(:,is_oce,k)
+              ENDIF 
+            END DO
+          ELSE 
+            DO i = 1, klon
+              pctsrf_t(i,is_ter,k) = pctsrf(i,is_ter)
+              IF (NINT(pctsrf(i,is_ter)).EQ.1 ) THEN
+                pctsrf_t(i,is_sic,k) = 0.
+                pctsrf_t(i,is_oce,k) = 0.                  
+                IF(phy_ice(i) .GE. 1.e-5) THEN
+                  pctsrf_t(i,is_lic,k) = phy_ice(i)
+                  pctsrf_t(i,is_ter,k) = pctsrf_t(i,is_ter,k) 
+     .                                   - pctsrf_t(i,is_lic,k)
+                ELSE
+                  pctsrf_t(i,is_lic,k) = 0.
+                ENDIF 
+              ELSE
+                pctsrf_t(i,is_lic,k) = 0. 
+                IF(phy_ice(i) .GE. 1.e-5) THEN 
+                  pctsrf_t(i,is_ter,k) = 0.
+                  pctsrf_t(i,is_sic,k) = phy_ice(i)
+                  pctsrf_t(i,is_oce,k) = 1. - pctsrf_t(i,is_sic,k)
+                ELSE
+                  pctsrf_t(i,is_sic,k) = 0.
+                  pctsrf_t(i,is_oce,k) = 1.                      
+                ENDIF 
+              ENDIF
+              verif = pctsrf_t(i,is_ter,k) +
+     .                pctsrf_t(i,is_oce,k) + 
+     .                pctsrf_t(i,is_sic,k) +
+     .                pctsrf_t(i,is_lic,k)
+              IF ( verif .LT. 1. - 1.e-5 .OR. 
+     $             verif .GT. 1 + 1.e-5) THEN  
+                WRITE(*,*) 'pb sous maille au point : i,k,verif '
      $                    , i,k,verif
-                  ENDIF 
-                END DO
-            ELSE 
-                DO i = 1, klon
-                  phy_nat(i,k) = phy_nat0(i)
-                  IF (NINT(phy_nat0(i)).EQ.1 ) THEN
-                      IF(phy_ice(i,k) .GE. 1.e-5) THEN 
-                          phy_icet(i,k) = phy_ice(i,k)
-                          phy_ice(i,k) = 0.
-                          phy_nat(i,k) = phy_nat(i,k) - phy_icet(i,k)
-                          phy_oce(i,k) = 0.
-                      ELSE
-                          phy_icet(i,k) = 0.
-                          phy_ice(i,k) = 0.
-                          phy_oce(i,k) = 0.                      
-                      ENDIF 
-                  ELSE 
-                      IF(phy_ice(i,k) .GE. 1.e-5) THEN 
-                          phy_icet(i,k) = 0.
-                          phy_nat(i,k) = 0.
-                          phy_oce(i,k) = 1. - phy_ice(i,k)
-                      ELSE
-                          phy_icet(i,k) = 0.
-                          phy_ice(i,k) = 0.
-                          phy_oce(i,k) = 1.                      
-                      ENDIF 
-                  ENDIF
-                  verif = phy_nat(i,k) + phy_icet(i,k)+ phy_ice(i,k) 
-     $                + phy_oce(i,k)
-                  IF ( verif .LT. 1. - 1.e-5 .OR. 
-     $                verif .GT. 1 + 1.e-5) THEN  
-                      WRITE(*,*) 'pb sous maille au point : i,k,verif '
-     $                    , i,k,verif
-                  ENDIF 
-                END DO
-            ENDIF 
+              ENDIF 
+            END DO
+          ENDIF 
         ELSE  
-            DO i = 1, klon
-              phy_nat(i,k) = phy_nat0(i)
-              IF ( (phy_ice(i,k) - 0.5).GE.1.e-5 ) THEN
-                  IF (NINT(phy_nat0(i)).EQ.0) THEN
-                      phy_nat(i,k) = 3.0
-                  ELSE
-                      phy_nat(i,k) = 2.0
-                  ENDIF
+          DO i = 1, klon
+            phy_nat(i,k) = phy_nat0(i)
+            IF ( (phy_ice(i) - 0.5).GE.1.e-5 ) THEN
+              IF (NINT(phy_nat0(i)).EQ.0) THEN
+                phy_nat(i,k) = 3.0
+              ELSE
+                phy_nat(i,k) = 2.0
               ENDIF
-            END DO
+            ENDIF
+          END DO
         ENDIF 
       ENDDO
@@ -892,11 +951,11 @@
       IF (newlmt ) THEN
           ierr = NF_PUT_VARA_DOUBLE (nid,id_FOCE,debut,epais
-     $        ,phy_oce(1,k))
+     $        ,pctsrf_t(1,is_oce,k))
           ierr = NF_PUT_VARA_DOUBLE (nid,id_FSIC,debut,epais
-     $        ,phy_ice(1,k))
+     $        ,pctsrf_t(1,is_sic,k))
           ierr = NF_PUT_VARA_DOUBLE (nid,id_FTER,debut,epais
-     $        ,phy_nat(1,k))
+     $        ,pctsrf_t(1,is_ter,k))
           ierr = NF_PUT_VARA_DOUBLE (nid,id_FLIC,debut,epais
-     $        ,phy_icet(1,k))
+     $        ,pctsrf_t(1,is_lic,k))
       ELSE 
           ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais
@@ -912,11 +971,11 @@
       IF (newlmt ) THEN
           ierr = NF_PUT_VARA_REAL (nid,id_FOCE,debut,epais
-     $        ,phy_oce(1,k))
+     $        ,pctsrf_t(1,is_oce,k))
           ierr = NF_PUT_VARA_REAL (nid,id_FSIC,debut,epais
-     $        ,phy_ice(1,k))
+     $        ,pctsrf_t(1,is_sic,k))
           ierr = NF_PUT_VARA_REAL (nid,id_FTER,debut,epais
-     $        ,phy_nat(1,k))
+     $        ,pctsrf_t(1,is_ter,k))
           ierr = NF_PUT_VARA_REAL (nid,id_FLIC,debut,epais
-     $        ,phy_icet(1,k))
+     $        ,pctsrf_t(1,is_lic,k))
       ELSE 
           ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais
