Changeset 5022

Timestamp:

Jul 5, 2024, 4:38:48 PM (2 weeks ago)

Author:

Sebastien Nguyen

Message:

include ISO keys in pbl_surface and associated routines in phylmd

Location:

LMDZ6/trunk/libf/phylmd

Files:

• fonte_neige_mod.F90 (modified) (17 diffs)
• ocean_forced_mod.F90 (modified) (16 diffs)
• pbl_surface_mod.F90 (modified) (68 diffs)
• phys_output_write_mod.F90 (modified) (4 diffs)
• surf_land_bucket_mod.F90 (modified) (8 diffs)
• surf_land_mod.F90 (modified) (13 diffs)
• surf_landice_mod.F90 (modified) (10 diffs)
• surf_ocean_mod.F90 (modified) (7 diffs)
• surf_seaice_mod.F90 (modified) (6 diffs)

LMDZ6/trunk/libf/phylmd/fonte_neige_mod.F90

@@ -36 +36 @@
   REAL, ALLOCATABLE, DIMENSION(:)             :: runofflic_global
   !$OMP THREADPRIVATE(runofflic_global)
+#ifdef ISO
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE  :: xtrun_off_ter
+  !$OMP THREADPRIVATE(xtrun_off_ter)
+  REAL, ALLOCATABLE, DIMENSION(:,:)           :: xtrun_off_lic
+  !$OMP THREADPRIVATE(xtrun_off_lic)
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE  :: xtrun_off_lic_0 
+  !$OMP THREADPRIVATE(xtrun_off_lic_0)
+  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE:: fxtfonte_global
+  !$OMP THREADPRIVATE(fxtfonte_global)
+  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE:: fxtcalving_global
+  !$OMP THREADPRIVATE(fxtcalving_global)
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE  :: xtrunofflic_global
+  !$OMP THREADPRIVATE(xtrunofflic_global)
+#endif
 
 CONTAINS
@@ -123 +137 @@
 
   END SUBROUTINE fonte_neige_init
+
+#ifdef ISO
+  SUBROUTINE fonte_neige_init_iso(xtrestart_runoff)
+
+! This subroutine allocates and initialize variables in the module. 
+! The variable run_off_lic_0 is initialized to the field read from
+! restart file. The other variables are initialized to zero.
+
+    USE infotrac_phy, ONLY: niso
+#ifdef ISOVERIF
+    USE isotopes_mod, ONLY: iso_eau,iso_HDO
+    USE isotopes_verif_mod
+#endif
+!
+!****************************************************************************************
+! Input argument
+    REAL, DIMENSION(niso,klon), INTENT(IN) :: xtrestart_runoff 
+
+! Local variables
+    INTEGER                           :: error
+    CHARACTER (len = 80)              :: abort_message 
+    CHARACTER (len = 20)              :: modname = 'fonte_neige_init'
+    INTEGER                           :: i
+
+
+!****************************************************************************************
+! Allocate run-off at landice and initilize with field read from restart
+!
+!****************************************************************************************
+
+    ALLOCATE(xtrun_off_lic_0(niso,klon), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation run_off_lic'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF    
+    
+    xtrun_off_lic_0(:,:) = xtrestart_runoff(:,:)        
+
+#ifdef ISOVERIF 
+      IF (iso_eau > 0) THEN   
+        CALL iso_verif_egalite_vect1D( &
+     &           xtrun_off_lic_0,run_off_lic_0,'fonte_neige 100', &
+     &           niso,klon)
+      ENDIF !IF (iso_eau > 0) THEN
+#endif        
+
+!****************************************************************************************
+! Allocate other variables and initilize to zero
+!
+!****************************************************************************************
+
+    ALLOCATE(xtrun_off_ter(niso,klon), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation xtrun_off_ter'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    xtrun_off_ter(:,:) = 0.
+    
+    ALLOCATE(xtrun_off_lic(niso,klon), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation xtrun_off_lic'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    xtrun_off_lic(:,:) = 0.
+
+    ALLOCATE(fxtfonte_global(niso,klon,nbsrf))
+    IF (error /= 0) THEN
+       abort_message='Pb allocation fxtfonte_global'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    fxtfonte_global(:,:,:) = 0.0
+
+    ALLOCATE(fxtcalving_global(niso,klon,nbsrf))
+    IF (error /= 0) THEN
+       abort_message='Pb allocation fxtcalving_global'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    fxtcalving_global(:,:,:) = 0.0
+
+    ALLOCATE(xtrunofflic_global(niso,klon))
+    IF (error /= 0) THEN
+       abort_message='Pb allocation xtrunofflic_global'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    xtrunofflic_global(:,:) = 0.0
+
+  END SUBROUTINE fonte_neige_init_iso
+#endif
+
 !
 !****************************************************************************************
@@ -128 +231 @@
   SUBROUTINE fonte_neige( knon, nisurf, knindex, dtime, &
        tsurf, precip_rain, precip_snow, &
-       snow, qsol, tsurf_new, evap)
-
-  USE indice_sol_mod
+       snow, qsol, tsurf_new, evap &
+#ifdef ISO    
+     & ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag   &
+     & ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag   &
+#endif
+     &   )
+
+    USE indice_sol_mod
+#ifdef ISO
+    USE infotrac_phy, ONLY: niso
+    !use isotopes_mod, ONLY: ridicule_snow,iso_eau,iso_HDO
+#ifdef ISOVERIF
+    USE isotopes_verif_mod
+#endif
+#endif
         
 ! Routine de traitement de la fonte de la neige dans le cas du traitement
@@ -172 +287 @@
     REAL, DIMENSION(klon), INTENT(INOUT) :: tsurf_new
     REAL, DIMENSION(klon), INTENT(INOUT) :: evap
+
+#ifdef ISO    
+        ! sortie de quelques diagnostiques
+    REAL, DIMENSION(klon), INTENT(OUT) :: fq_fonte_diag
+    REAL, DIMENSION(klon), INTENT(OUT) :: fqfonte_diag
+    REAL, DIMENSION(klon), INTENT(OUT) ::  snow_evap_diag 
+    REAL, DIMENSION(klon), INTENT(OUT) ::  fqcalving_diag  
+    REAL,                  INTENT(OUT) :: max_eau_sol_diag  
+    REAL, DIMENSION(klon), INTENT(OUT) ::  runoff_diag   
+    REAL, DIMENSION(klon), INTENT(OUT) :: run_off_lic_diag  
+    REAL,                  INTENT(OUT) :: coeff_rel_diag
+#endif
 
 ! Local variables
@@ -193 +320 @@
 
     LOGICAL               :: neige_fond
+
+#ifdef ISO
+        max_eau_sol_diag=max_eau_sol
+#endif
+
 
 !****************************************************************************************
@@ -231 +363 @@
     
     bil_eau_s(:) = (precip_rain(:) * dtime) - (evap(:) - snow_evap(:)) * dtime
+#ifdef ISO
+    snow_evap_diag(:) = snow_evap(:)
+    coeff_rel_diag    = coeff_rel
+#endif
+
 
 
@@ -254 +391 @@
           bil_eau_s(i) = bil_eau_s(i) + fq_fonte 
           tsurf_new(i) = tsurf_new(i) - fq_fonte * chasno  
+#ifdef ISO
+          fq_fonte_diag(i) = fq_fonte
+#endif
+
 
 !IM cf JLD OK     
@@ -273 +414 @@
        snow(i)=MIN(snow(i),snow_max)
     ENDDO
+#ifdef ISO
+    DO i = 1, knon
+       fqcalving_diag(i) = fqcalving(i)
+       fqfonte_diag(i)   = fqfonte(i)
+    ENDDO !DO i = 1, knon
+#endif
+
 
     IF (nisurf == is_ter) THEN
@@ -278 +426 @@
           qsol(i) = qsol(i) + bil_eau_s(i)
           run_off_ter(i) = run_off_ter(i) + MAX(qsol(i) - max_eau_sol, 0.0)
+#ifdef ISO
+          runoff_diag(i) = MAX(qsol(i) - max_eau_sol, 0.0)
+#endif
           qsol(i) = MIN(qsol(i), max_eau_sol) 
        ENDDO
@@ -290 +441 @@
        ENDDO
     ENDIF
+
+#ifdef ISO
+    DO i = 1, klon    
+      run_off_lic_diag(i) = run_off_lic(i)
+    ENDDO ! DO i = 1, knon    
+#endif
     
 !****************************************************************************************
@@ -312 +469 @@
 !****************************************************************************************
 !
-  SUBROUTINE fonte_neige_final(restart_runoff)
+  SUBROUTINE fonte_neige_final(restart_runoff &
+#ifdef ISO      
+     &                        ,xtrestart_runoff &
+#endif   
+     &                        )
 !
 ! This subroutine returns run_off_lic_0 for later writing to restart file.
 !
+#ifdef ISO
+    USE infotrac_phy, ONLY: niso
+#ifdef ISOVERIF
+    USE isotopes_mod, ONLY: iso_eau
+    USE isotopes_verif_mod
+#endif
+#endif
+!
 !****************************************************************************************
     REAL, DIMENSION(klon), INTENT(OUT) :: restart_runoff
+#ifdef ISO     
+    REAL, DIMENSION(niso,klon), INTENT(OUT) :: xtrestart_runoff
+#ifdef ISOVERIF
+    INTEGER :: i
+#endif  
+#endif
+
+
 
 !****************************************************************************************
 ! Set the output variables
     restart_runoff(:) = run_off_lic_0(:)
+#ifdef ISO
+    xtrestart_runoff(:,:) = xtrun_off_lic_0(:,:)
+#ifdef ISOVERIF 
+    IF (iso_eau > 0) THEN   
+      DO i=1,klon
+        IF (iso_verif_egalite_nostop(run_off_lic_0(i) &
+     &                              ,xtrun_off_lic_0(iso_eau,i) &
+     &                              ,'fonte_neige 413') &
+     &      == 1) then
+          WRITE(*,*) 'i=',i
+          STOP
+        ENDIF
+      ENDDO !DO i=1,klon
+    ENDIF !IF (iso_eau > 0) then 
+#endif    
+#endif
+
+
 
 ! Deallocation of all varaibles in the module
@@ -334 +529 @@
     IF (ALLOCATED(fqcalving_global)) DEALLOCATE(fqcalving_global)
     IF (ALLOCATED(runofflic_global)) DEALLOCATE(runofflic_global)
+#ifdef ISO
+    IF (ALLOCATED(xtrun_off_lic_0)) DEALLOCATE(xtrun_off_lic_0)
+    IF (ALLOCATED(xtrun_off_ter)) DEALLOCATE(xtrun_off_ter)
+    IF (ALLOCATED(xtrun_off_lic)) DEALLOCATE(xtrun_off_lic)
+    IF (ALLOCATED(fxtfonte_global)) DEALLOCATE(fxtfonte_global)
+    IF (ALLOCATED(fxtcalving_global)) DEALLOCATE(fxtcalving_global)
+    IF (ALLOCATED(xtrunofflic_global)) DEALLOCATE(xtrunofflic_global)
+#endif
+
 
   END SUBROUTINE fonte_neige_final
@@ -340 +544 @@
 !
   SUBROUTINE fonte_neige_get_vars(pctsrf, fqcalving_out, &
-       fqfonte_out, ffonte_out, run_off_lic_out)
+              fqfonte_out, ffonte_out, run_off_lic_out &
+#ifdef ISO     
+     &       ,fxtcalving_out, fxtfonte_out,xtrun_off_lic_out &
+#endif     
+     &       )
 
 
@@ -349 +557 @@
 !****************************************************************************************
 
-  USE indice_sol_mod
+    USE indice_sol_mod
+#ifdef ISO
+    USE infotrac_phy, ONLY: niso
+#endif
 
     REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
@@ -358 +569 @@
     REAL, DIMENSION(klon), INTENT(OUT)      :: run_off_lic_out
 
+#ifdef ISO
+    REAL, DIMENSION(niso,klon), INTENT(OUT) :: fxtcalving_out
+    REAL, DIMENSION(niso,klon), INTENT(OUT) :: fxtfonte_out
+    REAL, DIMENSION(niso,klon), INTENT(OUT) :: xtrun_off_lic_out
+    INTEGER   :: i,ixt
+#endif
+ 
     INTEGER   :: nisurf
 !****************************************************************************************
@@ -364 +582 @@
     fqfonte_out(:)   = 0.0
     fqcalving_out(:) = 0.0
+#ifdef ISO        
+    fxtfonte_out(:,:)   = 0.0
+    fxtcalving_out(:,:) = 0.0
+#endif
 
     DO nisurf = 1, nbsrf
@@ -373 +595 @@
     run_off_lic_out(:)=runofflic_global(:)
 
+#ifdef ISO
+    DO nisurf = 1, nbsrf
+      DO i=1,klon
+        DO ixt=1,niso
+          fxtfonte_out(ixt,i) = fxtfonte_out(ixt,i) + fxtfonte_global(ixt,i,nisurf)*pctsrf(i,nisurf)
+          fxtcalving_out(ixt,i) = fxtcalving_out(ixt,i) + fxtcalving_global(ixt,i,nisurf)*pctsrf(i,nisurf)
+        ENDDO !DO ixt=1,niso
+      ENDDO !DO i=1,klon
+    ENDDO !DO nisurf = 1, nbsrf
+    xtrun_off_lic_out(:,:) = xtrunofflic_global(:,:)
+#endif
+
   END SUBROUTINE fonte_neige_get_vars
 !
 !****************************************************************************************
 !
+!#ifdef ISO
+!  subroutine fonte_neige_export_xtrun_off_lic_0(knon,xtrun_off_lic_0_diag)
+!    use infotrac_phy, ONLY: niso
+!
+!    ! inputs
+!    INTEGER, INTENT(IN)                      :: knon
+!    real, INTENT(IN), DIMENSION(niso,klon)   :: xtrun_off_lic_0_diag
+!
+!    xtrun_off_lic_0(:,:)=xtrun_off_lic_0_diag(:,:)
+!
+!  end subroutine fonte_neige_export_xtrun_off_lic_0
+!#endif 
+
+#ifdef ISO
+  SUBROUTINE gestion_neige_besoin_varglob_fonte_neige(klon,knon, &
+     &           xtprecip_snow,xtprecip_rain, &
+     &           fxtfonte_neige,fxtcalving, &
+     &           knindex,nisurf,run_off_lic_diag,coeff_rel_diag)
+
+        ! dans cette routine, on a besoin des variables globales de
+        ! fonte_neige_mod. Il faut donc la mettre dans fonte_neige_mod
+        ! le reste de gestion_neige est dans isotopes_routines_mod car sinon pb
+        ! de dépendance circulaire.
+
+    USE infotrac_phy, ONLY: ntiso,niso
+    USE isotopes_mod, ONLY: iso_eau    
+    USE indice_sol_mod    
+#ifdef ISOVERIF
+    USE isotopes_verif_mod
+#endif
+    IMPLICIT NONE
+
+    ! inputs
+    INTEGER, INTENT(IN)                     :: klon,knon
+    REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtprecip_snow, xtprecip_rain
+    REAL, DIMENSION(niso,klon), INTENT(IN)  :: fxtfonte_neige,fxtcalving
+    INTEGER, INTENT(IN)                     :: nisurf
+    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
+    REAL, DIMENSION(klon), INTENT(IN)       :: run_off_lic_diag  
+    REAL, INTENT(IN)                        :: coeff_rel_diag  
+
+    ! locals
+    INTEGER :: i,ixt,j
+        
+#ifdef ISOVERIF
+    IF (nisurf == is_lic) THEN
+      IF (iso_eau > 0) THEN  
+        DO i = 1, knon
+           j = knindex(i)
+           CALL iso_verif_egalite(xtrun_off_lic_0(iso_eau,j), &
+     &             run_off_lic_0(j),'gestion_neige_besoin_varglob_fonte_neige 625')
+        ENDDO
+      ENDIF
+    ENDIF
+#endif
+
+! calcul de run_off_lic 
+
+    IF (nisurf == is_lic) THEN
+!         coeff_rel = dtime/(tau_calv * rday)
+
+      DO i = 1, knon
+        j = knindex(i)
+        DO ixt = 1, niso
+          xtrun_off_lic(ixt,i)   = (coeff_rel_diag *  fxtcalving(ixt,i)) &
+     &                            +(1. - coeff_rel_diag) * xtrun_off_lic_0(ixt,j)
+          xtrun_off_lic_0(ixt,j) = xtrun_off_lic(ixt,i)
+          xtrun_off_lic(ixt,i)   = xtrun_off_lic(ixt,i) + fxtfonte_neige(ixt,i) + xtprecip_rain(ixt,i)
+        ENDDO !DO ixt=1,niso
+#ifdef ISOVERIF
+          IF (iso_eau > 0) THEN             
+            IF (iso_verif_egalite_choix_nostop(xtrun_off_lic(iso_eau,i), &
+     &                  run_off_lic_diag(i),'gestion_neige_besoin_varglob_fonte_neige 1201a', &
+     &                  errmax,errmaxrel) == 1) THEN 
+               WRITE(*,*) 'i,j=',i,j   
+               WRITE(*,*) 'coeff_rel_diag=',coeff_rel_diag
+               STOP
+            ENDIF
+          ENDIF
+#endif
+      ENDDO
+    ENDIF !IF (nisurf == is_lic) THEN  
+
+! Save ffonte, fqfonte and fqcalving in global arrays for each 
+! sub-surface separately
+    DO i = 1, knon
+      DO ixt = 1, niso
+        fxtfonte_global(ixt,knindex(i),nisurf)   = fxtfonte_neige(ixt,i)
+        fxtcalving_global(ixt,knindex(i),nisurf) = fxtcalving(ixt,i)
+      ENDDO !do ixt=1,niso
+    ENDDO   
+
+    IF (nisurf == is_lic) THEN
+      DO i = 1, knon    
+        DO ixt = 1, niso   
+        xtrunofflic_global(ixt,knindex(i)) = xtrun_off_lic(ixt,i)
+        ENDDO ! DO ixt=1,niso   
+      ENDDO
+    ENDIF
+       
+  END SUBROUTINE gestion_neige_besoin_varglob_fonte_neige
+#endif
+
+
 END MODULE fonte_neige_mod

LMDZ6/trunk/libf/phylmd/ocean_forced_mod.F90

@@ -22 +22 @@
        radsol, snow, agesno, & 
        qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-       tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa)
+       tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa &
+#ifdef ISO
+       ,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, &
+       xtsnow,xtevap,h1 &  
+#endif            
+       )
 !
 ! This subroutine treats the "open ocean", all grid points that are not entierly covered
@@ -36 +41 @@
     USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
     use config_ocean_skin_m, only: activate_ocean_skin
+#ifdef ISO
+    USE infotrac_phy, ONLY: ntiso,niso
+    USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, calcul_iso_surf_sic_vectall    
+#ifdef ISOVERIF
+    USE isotopes_mod, ONLY: iso_eau,ridicule
+    !USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix
+    USE isotopes_verif_mod
+#endif
+#endif
 
     INCLUDE "YOMCST.h"
@@ -57 +71 @@
     real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)
 
+#ifdef ISO
+    REAL, DIMENSION(ntiso,klon), INTENT(IN)  :: xtprecip_rain, xtprecip_snow
+    REAL, DIMENSION(ntiso,klon), INTENT(IN)  :: xtspechum
+    REAL, DIMENSION(klon),       INTENT(IN)  :: rlat
+#endif
+
 ! In/Output arguments
 !****************************************************************************************
@@ -62 +82 @@
     REAL, DIMENSION(klon), INTENT(INOUT)     :: snow
     REAL, DIMENSION(klon), INTENT(INOUT)     :: agesno !? put to 0 in ocean
-  
+#ifdef ISO     
+    REAL, DIMENSION(niso,klon), INTENT(IN)   :: xtsnow
+    REAL, DIMENSION(niso,klon), INTENT(INOUT):: Roce
+#endif
+
 ! Output arguments
 !****************************************************************************************
@@ -72 +96 @@
     REAL, intent(out):: sens_prec_liq(:) ! (knon)
 
+#ifdef ISO     
+    REAL, DIMENSION(ntiso,klon), INTENT(OUT) :: xtevap ! isotopes in evaporation flux
+    REAL, DIMENSION(klon),       INTENT(OUT) :: h1 ! just a diagnostic, not useful for the simulation
+#endif
+
 ! Local variables
 !****************************************************************************************
@@ -80 +109 @@
     REAL, DIMENSION(klon)       :: u1_lay, v1_lay
     LOGICAL                     :: check=.FALSE.
-    REAL sens_prec_sol(knon) 
-    REAL, DIMENSION(klon) :: lat_prec_liq, lat_prec_sol    
+    REAL, DIMENSION(knon)       :: sens_prec_sol
+    REAL, DIMENSION(klon)       :: lat_prec_liq, lat_prec_sol    
+#ifdef ISO   
+    REAL, PARAMETER :: t_coup = 273.15      
+#endif
+
 
 !****************************************************************************************
@@ -87 +120 @@
 !****************************************************************************************
     IF (check) WRITE(*,*)' Entering ocean_forced_noice'
-    
+
+#ifdef ISO
+#ifdef ISOVERIF
+    DO i = 1, knon
+      IF (iso_eau > 0) THEN         
+        CALL iso_verif_egalite_choix(xtspechum(iso_eau,i), &
+     &                  spechum(i),'ocean_forced_mod 111', &
+     &                  errmax,errmaxrel)     
+        CALL iso_verif_egalite_choix(snow(i), &
+     &                  xtsnow(iso_eau,i),'ocean_forced_mod 117', &
+     &                  errmax,errmaxrel)
+      ENDIF !IF (iso_eau > 0) THEN
+    ENDDO !DO i=1,knon
+#endif      
+#endif 
+
 !****************************************************************************************
 ! 1)    
@@ -103 +151 @@
 
     else ! GCM
-      CALL limit_read_sst(knon,knindex,tsurf_lim)
+      CALL limit_read_sst(knon,knindex,tsurf_lim &
+#ifdef ISO
+     &     ,Roce,rlat &
+#endif     
+     &     )
     endif ! knon
 !sb--
@@ -161 +213 @@
          flux_u1, flux_v1)  
 
+#ifdef ISO     
+    CALL calcul_iso_surf_oce_vectall(klon, knon,t_coup, &
+     &    ps,tsurf_new,spechum,u1_lay, v1_lay, xtspechum, &
+     &    evap, Roce,xtevap,h1 &
+#ifdef ISOTRAC
+     &    ,knindex &
+#endif
+     &    )
+#endif         
+
+#ifdef ISO
+#ifdef ISOVERIF
+!          write(*,*) 'ocean_forced_mod 176: sortie de ocean_forced_noice'
+    IF (iso_eau > 0) THEN
+      DO i = 1, knon               
+        CALL iso_verif_egalite_choix(snow(i), &
+     &          xtsnow(iso_eau,i),'ocean_forced_mod 180', &
+     &          errmax,errmaxrel)
+      ENDDO ! DO j=1,knon
+    ENDIF !IF (iso_eau > 0) THEN
+#endif
+#endif   
+
   END SUBROUTINE ocean_forced_noice
 !
@@ -173 +248 @@
        radsol, snow, qsol, agesno, tsoil, &
        qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-       tsurf_new, dflux_s, dflux_l, rhoa)
+       tsurf_new, dflux_s, dflux_l, rhoa &
+#ifdef ISO
+       ,xtprecip_rain, xtprecip_snow, xtspechum,Roce, &
+       xtsnow, xtsol,xtevap,Rland_ice &  
+#endif            
+       )
 !
 ! This subroutine treats the ocean where there is ice. 
@@ -187 +267 @@
     USE indice_sol_mod
     USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
+#ifdef ISO
+    USE infotrac_phy, ONLY: niso, ntiso
+    USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, calcul_iso_surf_sic_vectall 
+#ifdef ISOVERIF
+    USE isotopes_mod, ONLY: iso_eau,ridicule
+    !USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix
+    USE isotopes_verif_mod
+#endif
+#endif
 
 !   INCLUDE "indicesol.h"
@@ -209 +298 @@
     REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness
     real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)
+#ifdef ISO
+    REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow
+    REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtspechum
+    REAL, DIMENSION(niso,klon),  INTENT(IN) :: Roce 
+    REAL, DIMENSION(niso,klon),  INTENT(IN) :: Rland_ice
+#endif
 
 ! In/Output arguments
@@ -216 +311 @@
     REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
     REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
+#ifdef ISO     
+    REAL, DIMENSION(niso,klon), INTENT(INOUT)     :: xtsnow
+    REAL, DIMENSION(niso,klon), INTENT(IN)        :: xtsol
+#endif
 
 ! Output arguments
@@ -226 +325 @@
     REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
+#ifdef ISO     
+    REAL, DIMENSION(ntiso,klon), INTENT(OUT)      :: xtevap
+#endif      
 
 ! Local variables
@@ -238 +340 @@
     REAL, DIMENSION(klon)       :: u0, v0
     REAL, DIMENSION(klon)       :: u1_lay, v1_lay
-    REAL sens_prec_liq(knon), sens_prec_sol (knon)   
+    REAL, DIMENSION(knon)       :: sens_prec_liq, sens_prec_sol
     REAL, DIMENSION(klon)       :: lat_prec_liq, lat_prec_sol    
 
+#ifdef ISO
+    REAL, PARAMETER :: t_coup = 273.15
+    REAL, DIMENSION(klon) :: fq_fonte_diag
+    REAL, DIMENSION(klon) :: fqfonte_diag
+    REAL, DIMENSION(klon) :: snow_evap_diag 
+    REAL, DIMENSION(klon) :: fqcalving_diag 
+    REAL, DIMENSION(klon) :: run_off_lic_diag
+    REAL :: coeff_rel_diag
+    REAL :: max_eau_sol_diag  
+    REAL, DIMENSION(klon) :: runoff_diag    
+    INTEGER IXT
+    REAL, DIMENSION(niso,klon) :: xtsnow_prec, xtsol_prec
+    REAL, DIMENSION(klon) :: snow_prec, qsol_prec  
+#endif
 
 !****************************************************************************************
@@ -307 +423 @@
 !
 !****************************************************************************************
+#ifdef ISO
+   ! verif
+#ifdef ISOVERIF
+    DO i = 1, knon
+      IF (iso_eau > 0) THEN
+        IF (snow(i) > ridicule) THEN
+          CALL iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
+   &              'interfsurf 964',errmax,errmaxrel)
+        ENDIF !IF ((snow(i) > ridicule)) THEN
+      ENDIF !IF (iso_eau > 0) THEN      
+    ENDDO !DO i=1,knon  
+#endif
+   ! end verif
+
+    DO i = 1, knon
+      snow_prec(i) = snow(i)
+      DO ixt = 1, niso
+      xtsnow_prec(ixt,i) = xtsnow(ixt,i)
+      ENDDO !DO ixt=1,niso
+      ! initialisation:
+      fq_fonte_diag(i) = 0.0
+      fqfonte_diag(i)  = 0.0
+      snow_evap_diag(i)= 0.0
+    ENDDO !DO i=1,knon
+#endif
+
+
     CALL fonte_neige( knon, is_sic, knindex, dtime, &
          tsurf_tmp, precip_rain, precip_snow, &
-         snow, qsol, tsurf_new, evap)
+         snow, qsol, tsurf_new, evap &
+#ifdef ISO    
+     &  ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag   &
+     &  ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag   &
+#endif
+     &   )
+
+
+#ifdef ISO
+! isotopes: tout est externalisé
+!#ifdef ISOVERIF
+!        write(*,*) 'ocean_forced_mod 377: call calcul_iso_surf_sic_vectall'
+!        write(*,*) 'klon,knon=',klon,knon
+!#endif
+    CALL calcul_iso_surf_sic_vectall(klon,knon, &
+     &          evap,snow_evap_diag,Tsurf_new,Roce,snow, &
+     &          fq_fonte_diag,fqfonte_diag,dtime,t_coup, &
+     &          precip_snow,xtprecip_snow,xtprecip_rain, snow_prec,xtsnow_prec, &
+     &          xtspechum,spechum,ps, &
+     &          xtevap,xtsnow,fqcalving_diag, &
+     &          knindex,is_sic,run_off_lic_diag,coeff_rel_diag,Rland_ice &
+     &   )
+#ifdef ISOVERIF
+        !write(*,*) 'ocean_forced_mod 391: sortie calcul_iso_surf_sic_vectall'
+    IF (iso_eau > 0) THEN
+      DO i = 1, knon  
+        CALL iso_verif_egalite_choix(snow(i), &
+     &           xtsnow(iso_eau,i),'ocean_forced_mod 396', &
+     &           errmax,errmaxrel)
+      ENDDO ! DO j=1,knon
+    ENDIF !IF (iso_eau > 0) then
+#endif 
+!#ifdef ISOVERIF
+#endif   
+!#ifdef ISO
     
 ! Calculation of albedo at snow (alb_neig) and update the age of snow (agesno)

LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90

@@ -33 +33 @@
                                   wx_pbl_check, wx_pbl_dts_check, wx_evappot
   use config_ocean_skin_m, only: activate_ocean_skin
+#ifdef ISO
+  USE infotrac_phy, ONLY: niso,ntraciso=>ntiso    
+#endif
 
   IMPLICIT NONE
@@ -49 +52 @@
   !$OMP THREADPRIVATE(ydTs0, ydqs0)
 
+#ifdef ISO
+  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: xtsnow   ! snow at surface
+  !$OMP THREADPRIVATE(xtsnow)
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: Rland_ice   ! snow at surface
+  !$OMP THREADPRIVATE(Rland_ice)  
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: Roce   ! snow at surface
+  !$OMP THREADPRIVATE(Roce)  
+#endif
+
   INTEGER, SAVE :: iflag_pbl_surface_t2m_bug
   !$OMP THREADPRIVATE(iflag_pbl_surface_t2m_bug)
@@ -178 +190 @@
 
   END SUBROUTINE pbl_surface_init
+
+#ifdef ISO
+  SUBROUTINE pbl_surface_init_iso(xtsnow_rst,Rland_ice_rst)
+
+! This routine should be called after the restart file has been read.
+! This routine initialize the restart variables and does some validation tests
+! for the index of the different surfaces and tests the choice of type of ocean.
+
+    USE indice_sol_mod
+    USE print_control_mod, ONLY: lunout
+#ifdef ISOVERIF
+    USE isotopes_mod, ONLY: iso_eau,ridicule
+    USE isotopes_verif_mod
+#endif
+    IMPLICIT NONE
+
+    INCLUDE "dimsoil.h"
+ 
+! Input variables
+!****************************************************************************************
+    REAL, DIMENSION(niso,klon, nbsrf), INTENT(IN)          :: xtsnow_rst
+    REAL, DIMENSION(niso,klon), INTENT(IN)          :: Rland_ice_rst
+  
+! Local variables
+!****************************************************************************************
+    INTEGER                       :: ierr
+    CHARACTER(len=80)             :: abort_message
+    CHARACTER(len = 20)           :: modname = 'pbl_surface_init'
+    integer i,ixt
+    
+!****************************************************************************************
+! Allocate and initialize module variables with fields read from restart file.
+!
+!****************************************************************************************    
+
+    ALLOCATE(xtsnow(niso,klon,nbsrf), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
+
+    ALLOCATE(Rland_ice(niso,klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
+
+    ALLOCATE(Roce(niso,klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
+
+    xtsnow(:,:,:)  = xtsnow_rst(:,:,:)
+    Rland_ice(:,:) = Rland_ice_rst(:,:)
+    Roce(:,:)      = 0.0
+
+#ifdef ISOVERIF 
+      IF (iso_eau >= 0) THEN
+         CALL iso_verif_egalite_vect2D( &
+     &           xtsnow,snow, &
+     &           'pbl_surface_mod 170',niso,klon,nbsrf)
+         DO i=1,klon  
+            IF (iso_eau >= 0) THEN  
+              CALL iso_verif_egalite(Rland_ice(iso_eau,i),1.0, &
+     &         'pbl_surf_mod 177')
+            ENDIF
+         ENDDO
+      ENDIF
+#endif
+
+  END SUBROUTINE pbl_surface_init_iso
+#endif
+
 !  
 !****************************************************************************************
@@ -241 +318 @@
 !FC
 !!!
-                        )
+#ifdef ISO
+     &   ,xtrain_f, xtsnow_f,xt, &
+     &   wake_dlxt,zxxtevap,xtevap, &
+     &   d_xt,d_xt_w,d_xt_x, &
+     &   xtsol,dflux_xt,zxxtsnow,zxfluxxt,flux_xt, &
+     &   h1_diag,runoff_diag,xtrunoff_diag &
+#endif      
+     &   )
 !****************************************************************************************
 ! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
@@ -316 +400 @@
     USE mod_grid_phy_lmdz,  ONLY : nbp_lon, nbp_lat, grid1dto2d_glo
     USE print_control_mod,  ONLY : prt_level,lunout
+#ifdef ISO
+  USE isotopes_mod, ONLY: Rdefault,iso_eau
+#ifdef ISOVERIF
+        USE isotopes_verif_mod
+#endif
+#ifdef ISOTRAC
+        USE isotrac_mod, only: index_iso
+#endif
+#endif
     USE ioipsl_getin_p_mod, ONLY : getin_p
     use phys_state_var_mod, only: ds_ns, dt_ns, delta_sst, delta_sal, dter, &
@@ -368 +461 @@
     REAL, DIMENSION(klon),        INTENT(IN)        :: gustiness ! gustiness
 
-    REAL, DIMENSION(klon),        INTENT(IN)        :: cldt    ! total cloud fraction
+    REAL, DIMENSION(klon),        INTENT(IN)        :: cldt    ! total cloud 
+
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(IN)        :: xt       ! water vapour (kg/kg)
+    REAL, DIMENSION(ntraciso,klon),        INTENT(IN)        :: xtrain_f  ! rain fall
+    REAL, DIMENSION(ntraciso,klon),        INTENT(IN)        :: xtsnow_f  ! snow fall
+#endif
 
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
@@ -381 +480 @@
     REAL, DIMENSION(klon),        INTENT(IN)        :: wake_dens
 !!!
-
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(IN)        :: wake_dlxt   
+#endif
 ! Input/Output variables
 !****************************************************************************************
@@ -450 +551 @@
     REAL, INTENT(OUT):: zcoefm(:, :, :) ! (klon, klev, nbsrf + 1)
     ! coef for turbulent diffusion of U and V (?), mean for each grid point
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon),        INTENT(OUT)       :: zxxtevap     ! water vapour flux at surface, positiv upwards
+    REAL, DIMENSION(ntraciso,klon, klev),  INTENT(OUT)       :: d_xt        ! change in water vapour
+    REAL, DIMENSION(klon),                 INTENT(OUT)       :: runoff_diag
+    REAL, DIMENSION(niso,klon),            INTENT(OUT)       :: xtrunoff_diag
+    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)       :: d_xt_w
+    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)       :: d_xt_x
+#endif
+
+
 
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
@@ -513 +624 @@
     REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_v     ! v wind tension (kg m/s)/(m**2 s) or Pascal
 !FC
-    REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT)  :: treedrg      ! tree drag (m)               
+    REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: treedrg  ! tree drag (m)               
+#ifdef ISO        
+    REAL, DIMENSION(niso,klon),   INTENT(OUT)       :: xtsol      ! water height in the soil (mm)
+    REAL, DIMENSION(ntraciso,klon, nbsrf)           :: xtevap     ! evaporation at surface
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: h1_diag    ! just diagnostic, not useful
+#endif
 
 
@@ -527 +643 @@
     REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_qbs   ! blowind snow vertical flux (kg/m**2
 
+#ifdef ISO   
+    REAL, DIMENSION(ntraciso,klon),              INTENT(OUT) :: dflux_xt    ! change of water vapour flux
+    REAL, DIMENSION(niso,klon),                  INTENT(OUT) :: zxxtsnow    ! snow at surface, mean for each grid point
+    REAL, DIMENSION(ntraciso,klon, klev),        INTENT(OUT) :: zxfluxxt    ! water vapour flux, mean for each grid point 
+    REAL, DIMENSION(ntraciso,klon, klev, nbsrf), INTENT(OUT) :: flux_xt     ! water vapour flux(latent flux) (kg/m**2/s)  
+#endif
 
 ! Martin
@@ -575 +697 @@
     REAL, DIMENSION(klon)              :: ysnow, yqsurf, yagesno, yqsol
     REAL, DIMENSION(klon)              :: yrain_f, ysnow_f, ybs_f
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon)     :: yxt1
+    REAL, DIMENSION(niso,klon)         :: yxtsnow, yxtsol   
+    REAL, DIMENSION(ntraciso,klon)     :: yxtrain_f, yxtsnow_f 
+    REAL, DIMENSION(klon)              :: yrunoff_diag
+    REAL, DIMENSION(niso,klon)         :: yxtrunoff_diag
+    REAL, DIMENSION(niso,klon)         :: yRland_ice    
+#endif
     REAL, DIMENSION(klon)              :: ysolsw, ysollw
     REAL, DIMENSION(klon)              :: yfder
@@ -583 +713 @@
     REAL, DIMENSION(klon)              :: y_flux_t1, y_flux_q1
     REAL, DIMENSION(klon)              :: y_dflux_t, y_dflux_q
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon)     ::  y_flux_xt1
+    REAL, DIMENSION(ntraciso,klon)     ::  y_dflux_xt
+#endif
     REAL, DIMENSION(klon)              :: y_flux_u1, y_flux_v1
     REAL, DIMENSION(klon)              :: y_flux_bs, y_flux0
@@ -610 +744 @@
     REAL, DIMENSION(klon)              :: AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0
     REAL, DIMENSION(klon)              :: AcoefH, AcoefQ, BcoefH, BcoefQ
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon)     :: AcoefXT, BcoefXT
+#endif
     REAL, DIMENSION(klon)              :: AcoefU, AcoefV, BcoefU, BcoefV
     REAL, DIMENSION(klon)              :: AcoefQBS, BcoefQBS
@@ -628 +765 @@
     REAL, DIMENSION(klon,klev)         :: yu, yv
     REAL, DIMENSION(klon,klev)         :: yt, yq, yqbs
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon)      :: yxtevap
+    REAL, DIMENSION(ntraciso,klon,klev) :: y_d_xt
+    REAL, DIMENSION(ntraciso,klon,klev) :: y_flux_xt
+    REAL, DIMENSION(ntraciso,klon,klev) :: yxt   
+#endif
     REAL, DIMENSION(klon,klev)         :: ypplay, ydelp
     REAL, DIMENSION(klon,klev)         :: delp
@@ -699 +842 @@
     REAL, DIMENSION(klon,klev)         :: Kcoef_hq_w, Kcoef_m_w, gama_h_w, gama_q_w
     REAL, DIMENSION(klon)              :: alf_1, alf_2, alf_1_x, alf_2_x, alf_1_w, alf_2_w
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon,klev)         :: yxt_x, yxt_w
+    REAL, DIMENSION(ntraciso,klon)              :: y_flux_xt1_x , y_flux_xt1_w   
+    REAL, DIMENSION(ntraciso,klon,klev)         :: y_flux_xt_x,y_d_xt_x,zxfluxxt_x
+    REAL, DIMENSION(ntraciso,klon,klev)         :: y_flux_xt_w,y_d_xt_w,zxfluxxt_w
+    REAL, DIMENSION(ntraciso,klon,klev,nbsrf)   :: flux_xt_x, flux_xt_w
+    REAL, DIMENSION(ntraciso,klon)              :: AcoefXT_x, BcoefXT_x
+    REAL, DIMENSION(ntraciso,klon)              :: AcoefXT_w, BcoefXT_w
+    REAL, DIMENSION(ntraciso,klon,klev)         :: CcoefXT, DcoefXT
+    REAL, DIMENSION(ntraciso,klon,klev)         :: CcoefXT_x, DcoefXT_x
+    REAL, DIMENSION(ntraciso,klon,klev)         :: CcoefXT_w, DcoefXT_w
+    REAL, DIMENSION(ntraciso,klon,klev)         :: gama_xt,gama_xt_x,gama_xt_w
+#endif
 !!!
 !!!jyg le 08/02/2012
@@ -891 +1047 @@
     REAL, DIMENSION(klon)              :: yrmu0
     ! Martin
-    REAL, DIMENSIOn(klon)              :: yri0
+    REAL, DIMENSION(klon)              :: yri0
 
     REAL, DIMENSION(klon):: ydelta_sst, ydelta_sal, yds_ns, ydt_ns, ydter, &
@@ -898 +1054 @@
     ! dt_ds, tkt, tks, taur, sss on ocean points
     REAL :: missing_val
+#ifdef ISO
+    REAL, DIMENSION(klon)       :: h1
+    INTEGER                     :: ixt
+!#ifdef ISOVERIF
+!    integer iso_verif_positif_nostop
+!#endif    
+#endif
+
 !****************************************************************************************
 ! End of declarations
@@ -925 +1089 @@
       iflag_split_ref = mod(iflag_pbl_split,10)
       iflag_split = iflag_split_ref
+
+#ifdef ISO      
+#ifdef ISOVERIF
+      DO i=1,klon
+        DO ixt=1,niso
+          CALL iso_verif_noNaN(xtsol(ixt,i),'pbl_surface 608')
+        ENDDO
+      ENDDO
+#endif
+#ifdef ISOVERIF
+      DO i=1,klon  
+        IF (iso_eau >= 0) THEN  
+          CALL iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, &
+     &         'pbl_surf_mod 585',errmax,errmaxrel)
+          CALL iso_verif_egalite_choix(xtsnow_f(iso_eau,i),snow_f(i), &
+     &         'pbl_surf_mod 594',errmax,errmaxrel)
+          IF (iso_verif_egalite_choix_nostop(xtsol(iso_eau,i),qsol(i), &
+     &         'pbl_surf_mod 596',errmax,errmaxrel) == 1) THEN
+                WRITE(*,*) 'i=',i
+                STOP
+          ENDIF
+          DO nsrf=1,nbsrf
+            CALL iso_verif_egalite_choix(xtsnow(iso_eau,i,nsrf),snow(i,nsrf), &
+     &         'pbl_surf_mod 598',errmax,errmaxrel)
+          ENDDO
+        ENDIF !IF (iso_eau >= 0) THEN   
+      ENDDO !DO i=1,knon  
+      DO k=1,klev
+        DO i=1,klon  
+          IF (iso_eau >= 0) THEN  
+            CALL iso_verif_egalite_choix(xt(iso_eau,i,k),q(i,k), &
+     &           'pbl_surf_mod 595',errmax,errmaxrel)
+          ENDIF !IF (iso_eau >= 0) THEN  
+        ENDDO !DO i=1,knon  
+      ENDDO !DO k=1,klev
+#endif
+#endif
+
 
 !****************************************************************************************
@@ -990 +1192 @@
       PRINT*,'WARNING : On impose qsol=',qsol0
       qsol(:)=qsol0
+#ifdef ISO
+      DO ixt=1,niso
+        xtsol(ixt,:)=qsol0*Rdefault(ixt)
+      ENDDO
+#ifdef ISOTRAC      
+      DO ixt=1+niso,ntraciso
+        xtsol(ixt,:)=qsol0*Rdefault(index_iso(ixt))
+      ENDDO
+#endif       
+#endif
     ENDIF
 !****************************************************************************************
@@ -1040 +1252 @@
  qsnow(:)=0. ; snowhgt(:)=0. ; to_ice(:)=0. ; sissnow(:)=0.
  runoff(:)=0.
+#ifdef ISO
+zxxtevap(:,:)=0.
+ d_xt(:,:,:)=0. 
+ d_xt_x(:,:,:)=0.
+ d_xt_w(:,:,:)=0.
+ flux_xt(:,:,:,:)=0. 
+! xtsnow(:,:,:)=0.! attention, xtsnow est l'équivalent de snow et non de qsnow
+ xtevap(:,:,:)=0.
+#endif
     IF (iflag_pbl<20.or.iflag_pbl>=30) THEN
        zcoefh(:,:,:) = 0.0
@@ -1129 +1350 @@
 !FC
 
+#ifdef ISO
+   yxtrain_f = 0.0 ; yxtsnow_f = 0.0
+   yxtsnow  = 0.0
+   yxt = 0.0
+   yxtsol = 0.0
+   flux_xt = 0.0
+   yRland_ice = 0.0
+!   d_xt = 0.0      
+   y_dflux_xt = 0.0  
+   dflux_xt=0.0 
+   y_d_xt_x=0.      ; y_d_xt_w=0.       
+#endif 
+
 ! >> PC
 !the yfields_out variable is defined in (klon,nbcf_out) even if it is used on
@@ -1155 +1389 @@
     fluxlat_x(:,:)=0. ;           fluxlat_w(:,:)=0. 
 !>jyg
+#ifdef ISO
+    flux_xt_x(:,:,:,:)=0. ;          flux_xt_w(:,:,:,:)=0.
+#endif
 !
 !jyg<
@@ -1454 +1691 @@
           yfluxbs(j)=0.0
           y_flux_bs(j) = 0.0
+!!!
+#ifdef ISO
+          DO ixt=1,ntraciso
+            yxtrain_f(ixt,j) = xtrain_f(ixt,i)
+            yxtsnow_f(ixt,j) = xtsnow_f(ixt,i)  
+          ENDDO
+          DO ixt=1,niso
+            yxtsnow(ixt,j)   = xtsnow(ixt,i,nsrf)
+          ENDDO    
+          !IF (nsrf == is_lic) THEN
+          DO ixt=1,niso
+            yRland_ice(ixt,j)= Rland_ice(ixt,i)  
+          ENDDO    
+          !endif !IF (nsrf == is_lic) THEN
+#ifdef ISOVERIF
+          IF (iso_eau >= 0) THEN
+              call iso_verif_egalite_choix(ysnow_f(j), &
+     &          yxtsnow_f(iso_eau,j),'pbl_surf_mod 862', &
+     &          errmax,errmaxrel)
+              call iso_verif_egalite_choix(ysnow(j), &
+     &          yxtsnow(iso_eau,j),'pbl_surf_mod 872', &
+     &          errmax,errmaxrel)
+          ENDIF
+#endif
+#ifdef ISOVERIF
+         DO ixt=1,ntraciso
+           call iso_verif_noNaN(yxtsnow_f(ixt,j),'pbl_surf_mod 921')
+         ENDDO
+#endif
+#endif
        ENDDO
 ! >> PC
@@ -1493 +1760 @@
              yq(j,k) = q(i,k)
              yqbs(j,k)=qbs(i,k)
+#ifdef ISO
+             DO ixt=1,ntraciso   
+               yxt(ixt,j,k) = xt(ixt,i,k)
+             ENDDO !DO ixt=1,ntraciso
+#endif
           ENDDO
         ENDDO
@@ -1510 +1782 @@
              yq_w(j,k) = q(i,k)+(1.-wake_s(i))*wake_dlq(i,k)
 !!!
+#ifdef ISO
+             DO ixt=1,ntraciso
+               yxt_x(ixt,j,k) = xt(ixt,i,k)-wake_s(i)*wake_dlxt(ixt,i,k)
+               yxt_w(ixt,j,k) = xt(ixt,i,k)+(1.-wake_s(i))*wake_dlxt(ixt,i,k)
+             ENDDO
+#endif
           ENDDO
         ENDDO
@@ -1565 +1843 @@
              i = ni(j)
              yqsol(j) = qsol(i)
+#ifdef ISO
+             DO ixt=1,niso
+               yxtsol(ixt,j) = xtsol(ixt,i)
+             ENDDO
+#endif
           ENDDO
        ENDIF
@@ -1859 +2142 @@
             Kcoef_hq, gama_q, gama_h, &
 !!!
-            AcoefH, AcoefQ, BcoefH, BcoefQ)
+            AcoefH, AcoefQ, BcoefH, BcoefQ &
+#ifdef ISO
+         &   ,yxt, CcoefXT, DcoefXT, gama_xt, AcoefXT, BcoefXT & 
+#endif               
+         &   )
        ELSE  !(iflag_split .eq.0)
         CALL climb_hq_down(knon, ycoefh_x, ypaprs, ypplay, &
@@ -1867 +2154 @@
             Kcoef_hq_x, gama_q_x, gama_h_x, &
 !!!
-            AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x)
+            AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x &
+#ifdef ISO
+         &   ,yxt_x, CcoefXT_x, DcoefXT_x, gama_xt_x, AcoefXT_x, BcoefXT_x & 
+#endif               
+         &   )
 !!!
        IF (prt_level >=10) THEN
@@ -1882 +2173 @@
             Kcoef_hq_w, gama_q_w, gama_h_w, &
 !!!
-            AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w)
+            AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w &
+#ifdef ISO
+         &   ,yxt_w, CcoefXT_w, DcoefXT_w, gama_xt_w, AcoefXT_w, BcoefXT_w & 
+#endif               
+         &   )
 !!!
        IF (prt_level >=10) THEN
@@ -1964 +2259 @@
          yt1(:) = yt(:,1)
          yq1(:) = yq(:,1)
+#ifdef ISO
+         yxt1(:,:) = yxt(:,:,1)
+#endif
+
        ELSE IF (iflag_split .ge. 1) THEN
+#ifdef ISO
+        call abort_gcm('pbl_surface_mod 2149','isos pas encore dans iflag_split=1',1)
+#endif
+
 !
 ! Cdragq computation
@@ -2126 +2429 @@
                yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
                y_flux_u1, y_flux_v1, &
-               yveget,ylai,yheight )
+               yveget,ylai,yheight   &
+#ifdef ISO
+         &      ,yxtrain_f, yxtsnow_f,yxt1, &
+         &      yxtsnow,yxtsol,yxtevap,h1, &
+         &      yrunoff_diag,yxtrunoff_diag,yRland_ice &
+#endif               
+         &      )
  
 !FC quid qd yveget ylai yheight ne sont pas definit
@@ -2156 +2465 @@
           ENDDO
       ENDIF
-     
+
+#ifdef ISOVERIF
+        DO j=1,knon
+          DO ixt=1,ntraciso
+            CALL iso_verif_noNaN(yxtevap(ixt,j), &
+         &      'pbl_surface 1056a: apres surf_land')
+          ENDDO
+          DO ixt=1,niso
+            CALL iso_verif_noNaN(yxtsol(ixt,j), &
+         &      'pbl_surface 1056b: apres surf_land')
+          ENDDO
+        ENDDO
+#endif
+#ifdef ISOVERIF
+!        write(*,*) 'pbl_surface_mod 1038: sortie surf_land'
+        DO j=1,knon
+          IF (iso_eau >= 0) THEN     
+                 CALL iso_verif_egalite(yxtsnow(iso_eau,j), &
+     &                                  ysnow(j),'pbl_surf_mod 1043')
+          ENDIF !if (iso_eau.gt.0) then
+        ENDDO !DO i=1,klon
+#endif
+    
        CASE(is_lic)
           ! Martin
@@ -2177 +2508 @@
                   ysnowhgt, yqsnow, ytoice, ysissnow, &
                   yalb3_new, yrunoff, &
-                  y_flux_u1, y_flux_v1)
+                  y_flux_u1, y_flux_v1 &
+#ifdef ISO
+                  &    ,yxtrain_f, yxtsnow_f,yxt1,yRland_ice &
+                  &    ,yxtsnow,yxtsol,yxtevap &
+#endif              
+                  &    )
              
              !jyg<
@@ -2199 +2535 @@
                 ENDDO
              ENDIF
-             
+
+#ifdef ISOVERIF
+             DO j=1,knon
+               DO ixt=1,ntraciso
+                 CALL iso_verif_noNaN(yxtevap(ixt,j), &
+                        &             'pbl_surface 1095a: apres surf_landice')
+               ENDDO
+                do ixt=1,niso
+                   call iso_verif_noNaN(yxtsol(ixt,j), &
+                        &      'pbl_surface 1095b: apres surf_landice')
+                enddo
+             enddo
+#endif
+#ifdef ISOVERIF
+             !write(*,*) 'pbl_surface_mod 1060: sortie surf_landice'
+             do j=1,knon
+               IF (iso_eau >= 0) THEN     
+                 CALL iso_verif_egalite(yxtsnow(iso_eau,j), &
+                        &               ysnow(j),'pbl_surf_mod 1064')
+               ENDIF !if (iso_eau >= 0) THEN
+             ENDDO !DO i=1,klon
+#endif
+            
           END IF
           
@@ -2216 +2574 @@
                y_flux_u1, y_flux_v1, ydelta_sst(:knon), ydelta_sal(:knon), &
                yds_ns(:knon), ydt_ns(:knon), ydter(:knon), ydser(:knon), &
-               ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss)
+               ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss &
+#ifdef ISO
+         &      ,yxtrain_f, yxtsnow_f,yxt1,Roce, &
+         &      yxtsnow,yxtevap,h1 &
+#endif               
+         &      )
       IF (prt_level >=10) THEN
           print *,'arg de surf_ocean: ycdragh ',ycdragh(1:knon)
@@ -2257 +2620 @@
 !albedo SB <<<
                ytsurf_new, y_dflux_t, y_dflux_q, &
-               y_flux_u1, y_flux_v1)
+               y_flux_u1, y_flux_v1 &
+#ifdef ISO
+         &      ,yxtrain_f, yxtsnow_f,yxt1,Roce, &
+         &      yxtsnow,yxtsol,yxtevap,Rland_ice &
+#endif               
+         &      )
           
 ! Special DICE MPL 05082013 puis BOMEX MPL 20150410
@@ -2265 +2633 @@
           y_flux_v1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yv(j,1)*ypplay(j,1)/RD/yt(j,1)
           ENDDO
-      ENDIF
+       ENDIF
+
+#ifdef ISOVERIF
+        DO j=1,knon
+          DO ixt=1,ntraciso
+            CALL iso_verif_noNaN(yxtevap(ixt,j), &
+         &                       'pbl_surface 1165a: apres surf_seaice')
+          ENDDO
+          DO ixt=1,niso
+            CALL iso_verif_noNaN(yxtsol(ixt,j), &
+         &      'pbl_surface 1165b: apres surf_seaice')
+          ENDDO
+        ENDDO
+#endif
+#ifdef ISOVERIF
+        !write(*,*) 'pbl_surface_mod 1077: sortie surf_seaice'
+        DO j=1,knon
+          IF (iso_eau >= 0) THEN     
+                 CALL iso_verif_egalite(yxtsnow(iso_eau,j), &
+     &                                  ysnow(j),'pbl_surf_mod 1106')
+          ENDIF !IF (iso_eau >= 0) THEN
+        ENDDO !DO i=1,klon
+#endif
 
        CASE DEFAULT
@@ -2335 +2725 @@
           y_flux_t1(j) =  yfluxsens(j)
           y_flux_q1(j) = -yevap(j)
+#ifdef ISO
+          y_flux_xt1(:,:) = -yxtevap(:,:)
+#endif
           ENDDO
         ENDIF ! (ok_flux_surf)
@@ -2350 +2743 @@
 
        IF (iflag_split .GE. 1) THEN
+#ifdef ISO
+        call abort_gcm('pbl_surface_mod 2607','isos pas encore dans iflag_split=1',1)
+#endif
+!
 !
          IF (nsrf .ne. is_oce) THEN
@@ -2567 +2964 @@
             Kcoef_hq, gama_q, gama_h, &
 !!!
-            y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))    
+            y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:) &
+#ifdef ISO
+        &    ,yxt,y_flux_xt1 &
+        &    ,AcoefXT,BcoefXT,CcoefXT,DcoefXT,gama_xt &
+        &    ,y_flux_xt(:,:,:),y_d_xt(:,:,:) &
+#endif
+        &    )    
        ELSE  !(iflag_split .eq.0)
         CALL climb_hq_up(knon, dtime, yt_x, yq_x, &
@@ -2576 +2979 @@
             Kcoef_hq_x, gama_q_x, gama_h_x, &
 !!!
-            y_flux_q_x(:,:), y_flux_t_x(:,:), y_d_q_x(:,:), y_d_t_x(:,:))    
+            y_flux_q_x(:,:), y_flux_t_x(:,:), y_d_q_x(:,:), y_d_t_x(:,:) &
+#ifdef ISO
+        &    ,yxt_x,y_flux_xt1_x &
+        &    ,AcoefXT_x,BcoefXT_x,CcoefXT_x,DcoefXT_x,gama_xt_x &
+        &    ,y_flux_xt_x(:,:,:),y_d_xt_x(:,:,:) &
+#endif
+        &    )    
 !
        CALL climb_hq_up(knon, dtime, yt_w, yq_w, &
@@ -2585 +2994 @@
             Kcoef_hq_w, gama_q_w, gama_h_w, &
 !!!
-            y_flux_q_w(:,:), y_flux_t_w(:,:), y_d_q_w(:,:), y_d_t_w(:,:))    
+            y_flux_q_w(:,:), y_flux_t_w(:,:), y_d_q_w(:,:), y_d_t_w(:,:) &
+#ifdef ISO
+        &    ,yxt_w,y_flux_xt1_w &
+        &    ,AcoefXT_w,BcoefXT_w,CcoefXT_w,DcoefXT_w,gama_xt_w &
+        &    ,y_flux_xt_w(:,:,:),y_d_xt_w(:,:,:) &
+#endif
+        &    )    
 !!!
        ENDIF  ! (iflag_split .eq.0)
@@ -2703 +3118 @@
              flux_u(i,k,nsrf) = y_flux_u(j,k)
              flux_v(i,k,nsrf) = y_flux_v(j,k)
+
+#ifdef ISO
+             DO ixt=1,ntraciso
+                y_d_xt(ixt,j,k)  = y_d_xt(ixt,j,k) * ypct(j)
+                flux_xt(ixt,i,k,nsrf) = y_flux_xt(ixt,j,k)
+             ENDDO ! DO ixt=1,ntraciso
+             h1_diag(i)=h1(j)
+#endif
+
            ENDDO
         ENDDO
+
+#ifdef ISO
+#ifdef ISOVERIF
+        if (iso_eau.gt.0) then
+         call iso_verif_egalite_vect2D( &
+                y_d_xt,y_d_q, &
+                'pbl_surface_mod 2600',ntraciso,klon,klev)
+        endif       
+#endif
+#endif
 
        ELSE  !(iflag_split .eq.0)
@@ -2722 +3156 @@
             flux_u_x(i,k,nsrf) = y_flux_u_x(j,k)
             flux_v_x(i,k,nsrf) = y_flux_v_x(j,k)
+
+#ifdef ISO
+            DO ixt=1,ntraciso
+              y_d_xt_x(ixt,j,k)  = y_d_xt_x(ixt,j,k) * ypct(j)
+              flux_xt_x(ixt,i,k,nsrf) = y_flux_xt_x(ixt,j,k)
+            ENDDO ! DO ixt=1,ntraciso
+#endif
           ENDDO
         ENDDO
@@ -2739 +3180 @@
             flux_u_w(i,k,nsrf) = y_flux_u_w(j,k)
             flux_v_w(i,k,nsrf) = y_flux_v_w(j,k)
+
+#ifdef ISO
+            DO ixt=1,ntraciso
+              y_d_xt_w(ixt,j,k)  = y_d_xt_w(ixt,j,k) * ypct(j)
+              flux_xt_w(ixt,i,k,nsrf) = y_flux_xt_w(ixt,j,k)
+            ENDDO ! do ixt=1,ntraciso
+#endif
+
           ENDDO
         ENDDO
@@ -2750 +3199 @@
             flux_u(i,k,nsrf) = flux_u_x(i,k,nsrf)+ywake_s(j)*(flux_u_w(i,k,nsrf)-flux_u_x(i,k,nsrf))
             flux_v(i,k,nsrf) = flux_v_x(i,k,nsrf)+ywake_s(j)*(flux_v_w(i,k,nsrf)-flux_v_x(i,k,nsrf))
+#ifdef ISO
+            DO ixt=1,ntraciso
+              flux_xt(ixt,i,k,nsrf) = flux_xt_x(ixt,i,k,nsrf)+ywake_s(j)*(flux_xt_w(ixt,i,k,nsrf)-flux_xt_x(ixt,i,k,nsrf))
+            ENDDO ! do ixt=1,ntraciso
+#endif
           ENDDO
         ENDDO
@@ -2807 +3261 @@
           dflux_t(i) = dflux_t(i) + y_dflux_t(j)*ypct(j)
           dflux_q(i) = dflux_q(i) + y_dflux_q(j)*ypct(j)
+#ifdef ISO
+        DO ixt=1,niso
+          xtsnow(ixt,i,nsrf) = yxtsnow(ixt,j)  
+        ENDDO
+        DO ixt=1,ntraciso
+          xtevap(ixt,i,nsrf) = - flux_xt(ixt,i,1,nsrf)
+          dflux_xt(ixt,i) = dflux_xt(ixt,i) + y_dflux_xt(ixt,j)*ypct(j)
+        ENDDO  
+        IF (nsrf == is_lic) THEN
+          DO ixt=1,niso
+            Rland_ice(ixt,i) = yRland_ice(ixt,j)  
+          ENDDO
+        ENDIF !IF (nsrf == is_lic) THEN     
+#ifdef ISOVERIF
+        IF (iso_eau.gt.0) THEN  
+          call iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, &
+     &         'pbl_surf_mod 1230',errmax,errmaxrel)
+        ENDIF !if (iso_eau.gt.0) then
+#endif        
+#endif
        ENDDO
 
@@ -2911 +3385 @@
              i = ni(j)
              qsol(i) = yqsol(j)
+#ifdef ISO
+             runoff_diag(i)=yrunoff_diag(j)   
+             DO ixt=1,niso
+               xtsol(ixt,i) = yxtsol(ixt,j)
+               xtrunoff_diag(ixt,i)=yxtrunoff_diag(ixt,j)
+             ENDDO
+#endif
           ENDDO
        ENDIF
@@ -2923 +3404 @@
           ENDDO
        ENDDO
-       
+
+#ifdef ISO
+#ifdef ISOVERIF
+       !write(*,*) 'pbl_surface 2858'
+       DO i = 1, klon
+         DO ixt=1,niso
+           call iso_verif_noNaN(xtsol(ixt,i),'pbl_surface 1405')
+         ENDDO
+       ENDDO
+#endif
+#ifdef ISOVERIF
+     IF (iso_eau.gt.0) THEN
+        call iso_verif_egalite_vect2D( &
+                y_d_xt,y_d_q, &
+                'pbl_surface_mod 1261',ntraciso,klon,klev)
+     ENDIF !if (iso_eau.gt.0) then
+#endif
+#endif
 !!! jyg le 07/02/2012
        IF (iflag_split .ge.1) THEN
@@ -2942 +3440 @@
            d_u_w(i,k) = d_u_w(i,k) + y_d_u_w(j,k)
            d_v_w(i,k) = d_v_w(i,k) + y_d_v_w(j,k)
+#ifdef ISO
+           DO ixt=1,ntraciso
+             d_xt_x(ixt,i,k) = d_xt_x(ixt,i,k) + y_d_xt_x(ixt,j,k)
+             d_xt_w(ixt,i,k) = d_xt_w(ixt,i,k) + y_d_xt_w(ixt,j,k)
+           ENDDO ! DO ixt=1,ntraciso
+#endif
+
 !
 !!           d_wake_dlt(i,k) = d_wake_dlt(i,k) + y_d_t_w(i,k)-y_d_t_x(i,k)
@@ -2957 +3462 @@
              d_t(i,k) = d_t(i,k) + y_d_t(j,k)
              d_q(i,k) = d_q(i,k) + y_d_q(j,k)
+#ifdef ISO
+             DO ixt=1,ntraciso
+               d_xt(ixt,i,k) = d_xt(ixt,i,k) + y_d_xt(ixt,j,k)
+             ENDDO !DO ixt=1,ntraciso
+#endif
              d_u(i,k) = d_u(i,k) + y_d_u(j,k)
              d_v(i,k) = d_v(i,k) + y_d_v(j,k)
@@ -2971 +3481 @@
          ENDDO
         ENDIF
+
+#ifdef ISO
+#ifdef ISOVERIF
+!        write(*,*) 'd_q,d_xt(iso_eau,554,19)=',d_q(554,19),d_xt(iso_eau,554,19)
+!        write(*,*) 'pbl_surface 2929: d_q,d_xt(iso_eau,2,1)=',d_q(2,1),d_xt(iso_eau,2,1)
+!        write(*,*) 'y_d_q,y_d_xt(iso_eau,2,1)=',y_d_q(2,1),y_d_xt(iso_eau,2,1)
+!        write(*,*) 'iso_eau.gt.0=',iso_eau.gt.0
+        call iso_verif_noNaN_vect2D( &
+     &           d_xt, &
+     &           'pbl_surface 1385',ntraciso,klon,klev)  
+     IF (iso_eau >= 0) THEN
+        call iso_verif_egalite_vect2D( &
+                y_d_xt,y_d_q, &
+                'pbl_surface_mod 2945',ntraciso,klon,klev)
+        call iso_verif_egalite_vect2D( &
+                d_xt,d_q, &
+                'pbl_surface_mod 1276',ntraciso,klon,klev)
+     ENDIF !IF (iso_eau >= 0) THEN
+#endif
+#endif
 
 !      print*,'Dans pbl OK4'
@@ -3358 +3888 @@
    iflag_split=iflag_split_ref
 
+#ifdef ISO
+#ifdef ISOVERIF
+!        write(*,*) 'pbl_surface tmp 3249: d_q,d_xt(iso_eau,2,1)=',d_q(2,1),d_xt(iso_eau,2,1)
+    IF (iso_eau >= 0) THEN
+        call iso_verif_egalite_vect2D( &
+                d_xt,d_q, &
+                'pbl_surface_mod 1276',ntraciso,klon,klev)
+    ENDIF !IF (iso_eau >= 0) THEN
+#endif
+#endif
+
 !****************************************************************************************
 ! 16) Calculate the mean value over all sub-surfaces for some variables
@@ -3379 +3920 @@
     zxfluxt_w(:,:) = 0.0 ; zxfluxq_w(:,:) = 0.0
     zxfluxu_w(:,:) = 0.0 ; zxfluxv_w(:,:) = 0.0
+#ifdef ISO
+      zxfluxxt(:,:,:) = 0.0 
+      zxfluxxt_x(:,:,:) = 0.0
+      zxfluxxt_w(:,:,:) = 0.0
+#endif
+
 
 !!! jyg le 07/02/2012
@@ -3397 +3944 @@
               zxfluxu_w(i,k) = zxfluxu_w(i,k) + flux_u_w(i,k,nsrf) * pctsrf(i,nsrf)
               zxfluxv_w(i,k) = zxfluxv_w(i,k) + flux_v_w(i,k,nsrf) * pctsrf(i,nsrf)
+#ifdef ISO
+              DO ixt=1,ntraciso
+                zxfluxxt_x(ixt,i,k) = zxfluxxt_x(ixt,i,k) + flux_xt_x(ixt,i,k,nsrf) * pctsrf(i,nsrf)
+                zxfluxxt_w(ixt,i,k) = zxfluxxt_w(ixt,i,k) + flux_xt_w(ixt,i,k,nsrf) * pctsrf(i,nsrf)
+              ENDDO ! DO ixt=1,ntraciso
+#endif
             ENDDO
           ENDDO
@@ -3416 +3969 @@
              zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
              zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
+#ifdef ISO
+             DO ixt=1,niso
+               zxfluxxt(ixt,i,k) = zxfluxxt(ixt,i,k) + flux_xt(ixt,i,k,nsrf) * pctsrf(i,nsrf)
+             ENDDO ! DO ixt=1,niso
+#endif
           ENDDO
        ENDDO
@@ -3440 +3998 @@
        END DO
     endif
+
+#ifdef ISO
+    DO i = 1, klon
+      DO ixt=1,ntraciso
+        zxxtevap(ixt,i)     = - zxfluxxt(ixt,i,1)
+      ENDDO
+    ENDDO
+#endif
 
 !!!
@@ -3615 +4181 @@
     zxqsurf(:) = 0.0
     zxsnow(:)  = 0.0
+#ifdef ISO
+    zxxtsnow(:,:)  = 0.0
+#endif
+
     DO nsrf = 1, nbsrf
        DO i = 1, klon
           zxqsurf(i) = zxqsurf(i) + MAX(qsurf(i,nsrf),0.0) * pctsrf(i,nsrf)
           zxsnow(i)  = zxsnow(i)  + snow(i,nsrf)  * pctsrf(i,nsrf)
+#ifdef ISO
+          DO ixt=1,niso
+            zxxtsnow(ixt,i)  = zxxtsnow(ixt,i)  + xtsnow(ixt,i,nsrf)  * pctsrf(i,nsrf)
+          ENDDO ! DO ixt=1,niso
+#endif
        ENDDO
     ENDDO
@@ -3630 +4205 @@
 !****************************************************************************************
 !
-  SUBROUTINE pbl_surface_final(fder_rst, snow_rst, qsurf_rst, ftsoil_rst)
+  SUBROUTINE pbl_surface_final(fder_rst, snow_rst, qsurf_rst, ftsoil_rst &
+#ifdef ISO
+       ,xtsnow_rst,Rland_ice_rst &
+#endif       
+       )
 
     USE indice_sol_mod
+#ifdef ISO
+#ifdef ISOVERIF
+    USE isotopes_mod, ONLY: iso_eau,ridicule
+    USE isotopes_verif_mod, ONLY: errmax,errmaxrel
+#endif    
+#endif
 
     INCLUDE "dimsoil.h"
@@ -3642 +4227 @@
     REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: qsurf_rst
     REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
+#ifdef ISO
+    REAL, DIMENSION(niso,klon, nbsrf), INTENT(OUT)     :: xtsnow_rst
+    REAL, DIMENSION(niso,klon), INTENT(OUT)            :: Rland_ice_rst
+#endif
 
  
@@ -3652 +4241 @@
     qsurf_rst(:,:)    = qsurf(:,:)
     ftsoil_rst(:,:,:) = ftsoil(:,:,:)
+#ifdef ISO
+    xtsnow_rst(:,:,:)  = xtsnow(:,:,:) 
+    Rland_ice_rst(:,:) = Rland_ice(:,:)
+#endif
 
 !****************************************************************************************
@@ -3664 +4257 @@
     IF (ALLOCATED(ydTs0)) DEALLOCATE(ydTs0)
     IF (ALLOCATED(ydqs0)) DEALLOCATE(ydqs0)
+#ifdef ISO
+    IF (ALLOCATED(xtsnow)) DEALLOCATE(xtsnow)
+    IF (ALLOCATED(Rland_ice)) DEALLOCATE(Rland_ice)
+    IF (ALLOCATED(Roce)) DEALLOCATE(Roce)
+#endif
 
 !jyg<
@@ -3682 +4280 @@
   SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, &
        evap, z0m, z0h, agesno,                                  &
-       tsurf,alb_dir,alb_dif, ustar, u10m, v10m, tke)  
+       tsurf,alb_dir,alb_dif, ustar, u10m, v10m, tke &
+#ifdef ISO
+      ,xtevap  &
+#endif
+&      )  
     !albedo SB <<<
     ! Give default values where new fraction has appread
@@ -3711 +4313 @@
     REAL, DIMENSION(klon,nbsrf+1), INTENT(INOUT)        :: z0m,z0h
     REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(INOUT) :: tke
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon,nbsrf), INTENT(INOUT)        :: xtevap
+#endif
 
 ! Local variables
@@ -3718 +4323 @@
     CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
     INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
+#ifdef ISO
+    INTEGER           :: ixt
+#endif
 !
 ! All at once !! 
@@ -3763 +4371 @@
                 u10m(i,nsrf)  = u10m(i,nsrf_comp1)
                 v10m(i,nsrf)  = v10m(i,nsrf_comp1)
+#ifdef ISO
+                DO ixt=1,ntraciso
+                  xtevap(ixt,i,nsrf) = xtevap(ixt,i,nsrf_comp1)       
+                ENDDO       
+#endif
                 IF (iflag_pbl > 1) THEN
                  tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
@@ -3818 +4431 @@
                 u10m(i,nsrf)  = u10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
                 v10m(i,nsrf)  = v10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
+#ifdef ISO
+                DO ixt=1,ntraciso
+                  xtevap(ixt,i,nsrf) = xtevap(ixt,i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) &
+                                     + xtevap(ixt,i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3) 
+                ENDDO       
+#endif
                 IF (iflag_pbl > 1) THEN
                  tke(i,:,nsrf) = tke(i,:,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
@@ -3830 +4449 @@
              agesno(i,nsrf)   = 0.
              ftsoil(i,:,nsrf) = tsurf(i,nsrf)
+#ifdef ISO            
+             xtsnow(:,i,nsrf) = 0.
+#endif
           ELSE
              pfois(nsrf) = pfois(nsrf)+ 1

LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90

@@ -207 +207 @@
 ! Isotopes
          o_xtprecip,o_xtplul,o_xtpluc,o_xtovap,o_xtoliq,o_xtcond, &
+         o_runoff_diag, o_xtrunoff_diag, &
          o_xtevap,o_dxtdyn,o_dxtldyn,o_dxtcon,o_dxtlsc,o_dxteva, &
          o_dxtajs,o_dxtvdf,o_dxtthe, o_dxtch4, &
@@ -392 +393 @@
         d_xt_ajs, d_xt_ajsb, &
         d_xt_prod_nucl,d_xt_cosmo,d_xt_decroiss, &
+        runoff_diag, xtrunoff_diag, &
 #endif
          ep, epmax_diag, &  ! epmax_cape
@@ -2799 +2801 @@
 #ifdef ISO
     !write(*,*) 'tmp phys_output_write: ntiso=',ntiso
+    !! runoff land bucket - ajout S. Nguyen 25 avril 2024
+    CALL histwrite_phy(o_runoff_diag, runoff_diag)
+
     do ixt=1,ntiso
         !write(*,*) 'ixt,o_xtovap(ixt)=',ixt,o_xtovap(ixt)
@@ -2812 +2817 @@
         CALL histwrite_phy(o_xtovap(ixt),  xt_seri(ixt,:,:))
         CALL histwrite_phy(o_xtoliq(ixt), xtl_seri(ixt,:,:))
+
+        !! runoff land bucket - ajout S. Nguyen 25 avril 2024
+        CALL histwrite_phy(o_xtrunoff_diag(ixt), xtrunoff_diag(ixt,:))
+
 
         DO nsrf = 1, nbsrf ! ajout Camille 8 mai 2023

LMDZ6/trunk/libf/phylmd/surf_land_bucket_mod.F90

@@ -16 +16 @@
        snow, qsol, agesno, tsoil, &
        qsurf, z0_new, alb1_new, alb2_new, evap, &
-       fluxsens, fluxlat, tsurf_new, dflux_s, dflux_l)
+       fluxsens, fluxlat, tsurf_new, dflux_s, dflux_l &
+#ifdef ISO
+       ,xtprecip_rain, xtprecip_snow,xtspechum, &
+       xtsnow, xtsol,xtevap,h1, &
+       runoff_diag,xtrunoff_diag,Rland_ice &
+#endif           
+            )
 
     USE limit_read_mod
@@ -28 +34 @@
     USE mod_phys_lmdz_para
     USE indice_sol_mod
+#ifdef ISO
+    use infotrac_phy, ONLY: ntiso,niso
+    USE isotopes_mod, ONLY: iso_eau, iso_HDO, iso_O18, iso_O17, &
+        ridicule_qsol
+    USE isotopes_routines_mod, ONLY: calcul_iso_surf_ter_vectall
+#ifdef ISOVERIF
+    USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_noNaN, &
+        iso_verif_aberrant_o17,iso_verif_egalite_choix,iso_verif_egalite
+#endif
+#endif
 !****************************************************************************************
 ! Bucket calculations for surface. 
@@ -52 +68 @@
     REAL, DIMENSION(klon), INTENT(IN)       :: rugoro
     REAL, DIMENSION(klon), INTENT(IN)       :: swnet, lwnet
+#ifdef ISO
+    REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow
+    REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtspechum    
+#endif
 
 ! In/Output variables
@@ -58 +78 @@
     REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
     REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
+#ifdef ISO
+    REAL, DIMENSION(niso,klon), INTENT(INOUT)       :: xtsnow,xtsol
+#endif
 
 ! Output variables
@@ -67 +90 @@
     REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
+#ifdef ISO
+    REAL, DIMENSION(ntiso,klon), INTENT(OUT) :: xtevap
+    REAL, DIMENSION(klon),       INTENT(OUT) :: h1
+    REAL, DIMENSION(niso,klon),  INTENT(OUT) :: xtrunoff_diag
+    REAL, DIMENSION(klon),       INTENT(OUT) :: runoff_diag
+    REAL, DIMENSION(niso,klon),  INTENT(IN)  :: Rland_ice
+#endif
 
 ! Local variables
@@ -78 +108 @@
     REAL, DIMENSION(klon) :: dummy_riverflow,dummy_coastalflow 
     INTEGER               :: i
-!
-!****************************************************************************************
-
+#ifdef ISO
+    INTEGER               :: ixt
+    REAL, DIMENSION(niso,klon) :: xtsnow_prec,xtsol_prec
+    REAL, DIMENSION(klon) :: snow_prec,qsol_prec
+    REAL, PARAMETER       :: t_coup = 273.15
+    REAL, DIMENSION(klon) :: fq_fonte_diag
+    REAL, DIMENSION(klon) :: fqfonte_diag
+    REAL, DIMENSION(klon) :: snow_evap_diag 
+    REAL, DIMENSION(klon) :: fqcalving_diag 
+    REAL                  :: max_eau_sol_diag  
+    REAL, DIMENSION(klon) :: run_off_lic_diag 
+    REAL :: coeff_rel_diag
+#endif 
+!
+!****************************************************************************************
+
+#ifdef ISO
+#ifdef ISOVERIF
+        !write(*,*) 'surf_land_bucket 152'
+        DO i=1,knon
+          IF (iso_eau > 0) THEN
+            CALL iso_verif_egalite_choix(precip_snow(i), &
+     &                                   xtprecip_snow(iso_eau,i),'surf_land_bucket 131', &
+     &                                   errmax,errmaxrel)
+            CALL iso_verif_egalite_choix(qsol(i), &
+     &                                   xtsol(iso_eau,i),'surf_land_bucket 134', &
+     &                                   errmax,errmaxrel)
+          ENDIF
+        ENDDO
+#endif 
+#ifdef ISOVERIF
+        DO i=1,knon
+          DO ixt=1,niso
+            CALL iso_verif_noNaN(xtsol(ixt,i),'surf_land_mod_bucket 142')
+          ENDDO !do ixt=1,niso
+        ENDDO !do i=1,knon
+        !write(*,*) 'surf_land_bucket 152'
+#endif
+#endif
 
 !
@@ -131 +197 @@
          tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
     
+#ifdef ISO
+   ! verif
+#ifdef ISOVERIF
+    !write(*,*) 'surf_land_bucket 211'
+    DO i=1,knon
+      IF (iso_eau > 0) THEN
+        CALL iso_verif_egalite_choix(xtsnow(iso_eau,i), &
+     &           snow(i),'surf_land_bucket 522', &
+     &           errmax,errmaxrel) 
+      ENDIF !IF (iso_eau > 0) then
+    ENDDO !DO i=1,knon 
+#endif
+   ! end verif
+#endif         
+#ifdef ISO
+    DO i=1,knon
+      snow_prec(i)=snow(i)
+      qsol_prec(i)=qsol(i)
+      DO ixt=1,niso
+        xtsnow_prec(ixt,i)=xtsnow(ixt,i)
+        xtsol_prec(ixt,i) =xtsol(ixt,i)
+      ENDDO !DO ixt=1,niso
+      ! initialisation:
+      fqfonte_diag(i)  =0.0
+      fq_fonte_diag(i) =0.0
+      snow_evap_diag(i)=0.0
+    ENDDO !DO i=1,knon 
+#ifdef ISOVERIF
+    ! write(*,*) 'surf_land_bucket 235'
+    DO i=1,knon  
+      IF (iso_eau > 0) THEN
+        CALL iso_verif_egalite(qsol_prec(i),xtsol_prec(iso_eau,i), &
+    &                              'surf_land_bucket 141')
+      ENDIF
+    ENDDO !DO i=1,knon
+#endif    
+#endif    
 !
 !* Calculate snow height, run_off, age of snow
@@ -136 +239 @@
     CALL fonte_neige( knon, is_ter, knindex, dtime, &
          tsurf, precip_rain, precip_snow, &
-         snow, qsol, tsurf_new, evap)
+         snow, qsol, tsurf_new, evap &
+#ifdef ISO    
+     & ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag   &
+     & ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag   &
+#endif
+     &   )
+
+#ifdef ISO
+#ifdef ISOVERIF
+        DO i=1,knon
+          DO ixt=1,niso
+            CALL iso_verif_noNaN(xtsol_prec(ixt,i),'surf_land_burcket 237')
+          ENDDO
+        ENDDO
+#endif
+#ifdef ISOVERIF
+        !write(*,*) 'surf_land_bucket 235'
+        DO i=1,knon
+          IF (iso_eau > 0) THEN
+            CALL iso_verif_egalite_choix(qsol_prec(i), &
+     &                                   xtsol_prec(iso_eau,i),'surf_land_bucket 628', &
+     &                                   errmax,errmaxrel)
+            CALL iso_verif_egalite_choix(precip_snow(i), &
+     &                                   xtprecip_snow(iso_eau,i),'surf_land_bucket 227', &
+     &                                   errmax,errmaxrel)
+             ! attention, dans fonte_neige, on modifie snow sans modifier
+             ! xtsnow
+             ! c'est fait plus tard dans gestion_neige
+!            write(*,*) 'surf_land_bucket 287: i=',i
+!            write(*,*) 'snow(i)=',snow(i)
+            CALL iso_verif_egalite_choix(xtsnow(iso_eau,i), &
+     &                                   snow_prec(i),'surf_land_bucket 245', &
+     &                                   errmax,errmaxrel) 
+          ENDIF  
+          IF ((iso_O17 > 0).AND.(iso_O18 > 0)) THEN
+              IF (qsol_prec(i) > ridicule_qsol) THEN
+                CALL iso_verif_aberrant_o17(xtsol_prec(iso_O17,i)/qsol_prec(i) &
+     &                                     ,xtsol_prec(iso_O18,i)/qsol_prec(i) &
+     &                                     ,'surf_land_bucket 642')
+              ENDIF !IF ((qsol_prec(i) > ridicule_qsol) &
+          ENDIF !IF ((iso_O17 > 0).AND.(iso_O18 > 0)) THEN
+        ENDDO  !DO i=1,knon 
+        !write(*,*) 'surf_land_mod 291'
+        !write(*,*) 'snow_evap_diag(1)=',snow_evap_diag(1)
+#endif          
+        CALL calcul_iso_surf_ter_vectall(klon,knon, &
+     &           evap,snow_evap_diag,snow, &
+     &           fq_fonte_diag,fqfonte_diag,dtime,precip_rain,xtprecip_rain, &
+     &           precip_snow,xtprecip_snow, snow_prec,xtsnow_prec, &
+     &           tsurf_new,xtspechum,pref,spechum,t_coup,u1_lay,v1_lay,p1lay, &
+     &           qsol,xtsol,qsol_prec,xtsol_prec, &
+     &           max_eau_sol_diag, &
+     &           xtevap,xtsnow,h1,runoff_diag,xtrunoff_diag,fqcalving_diag, &
+     &           knindex,is_ter,run_off_lic_diag,coeff_rel_diag,Rland_ice &
+     &   )
+!#ifdef ISOVERIF
+!        write(*,*) 'surf_land_bucket 303'
+!#endif
+#endif
+
 !
 !* Calculate the age of snow

LMDZ6/trunk/libf/phylmd/surf_land_mod.F90

@@ -20 +20 @@
        qsurf, tsurf_new, dflux_s, dflux_l, &
        flux_u1, flux_v1 , & 
-       veget,lai,height)
+       veget,lai,height &
+#ifdef ISO
+       ,xtprecip_rain, xtprecip_snow,xtspechum, &
+       xtsnow, xtsol,xtevap,h1, &
+       runoff_diag,xtrunoff_diag,Rland_ice &
+#endif               
+               )
 
     USE dimphy
@@ -59 +65 @@
     USE calcul_fluxs_mod
     USE indice_sol_mod
+#ifdef ISO
+    use infotrac_phy, ONLY: ntiso,niso
+    use isotopes_mod, ONLY: nudge_qsol, iso_eau
+#ifdef ISOVERIF
+    use isotopes_verif_mod
+#endif
+#endif
+
     USE print_control_mod, ONLY: lunout
 
@@ -92 +106 @@
                                                          ! corresponds to previous sollwdown
     REAL, DIMENSION(klon), INTENT(IN)       :: q2m, t2m
-
+#ifdef ISO
+    REAL, DIMENSION(ntiso,klon), INTENT(IN)       :: xtprecip_rain, xtprecip_snow
+    REAL, DIMENSION(ntiso,klon), INTENT(IN)       :: xtspechum
+#endif 
 ! In/Output variables
 !****************************************************************************************
@@ -98 +115 @@
     REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
     REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
+#ifdef ISO
+    REAL, DIMENSION(niso,klon), INTENT(INOUT)    :: xtsnow, xtsol
+#endif
 
 ! Output variables
@@ -116 +136 @@
     REAL, DIMENSION(klon,nvm_lmdz), INTENT(OUT) :: veget,lai
     REAL, DIMENSION(klon,nvm_lmdz), INTENT(OUT) :: height
+#ifdef ISO
+    REAL, DIMENSION(ntiso,klon), INTENT(OUT)      :: xtevap
+    REAL, DIMENSION(klon), INTENT(OUT)      :: h1
+    REAL, DIMENSION(niso,klon), INTENT(OUT)      :: xtrunoff_diag
+    REAL, DIMENSION(klon), INTENT(OUT)      :: runoff_diag
+    REAL, DIMENSION(niso,klon), INTENT(IN)        :: Rland_ice
+#endif
 
 ! Local variables
@@ -132 +159 @@
 !albedo SB <<<
 
-
+#ifdef ISO       
+      real, parameter :: t_coup = 273.15
+      real, dimension(klon) :: fqfonte_diag
+      real, dimension(klon) :: snow_evap_diag 
+      real, dimension(klon) :: fqcalving_diag 
+      integer :: ixt
+#endif 
 !****************************************************************************************
 !Total solid precip
@@ -142 +175 @@
 ENDIF
 !****************************************************************************************
+#ifdef ISO
+#ifdef ISOVERIF
+!        write(*,*) 'surf_land_mod 162'
+        do i=1,knon
+          if (iso_eau.gt.0) then
+            call iso_verif_egalite_choix(precip_snow(i), &
+     &          xtprecip_snow(iso_eau,i),'surf_land_mod 129', &
+     &          errmax,errmaxrel)
+            call iso_verif_egalite_choix(qsol(i), &
+     &          xtsol(iso_eau,i),'surf_land_mod 139', &
+     &          errmax,errmaxrel)
+          endif  
+        enddo
+#endif 
+#ifdef ISOVERIF
+!       write(*,*) 'surf_land 169: ok_veget=',ok_veget
+        do i=1,knon
+         do ixt=1,ntiso
+           call iso_verif_noNaN(xtprecip_snow(ixt,i),'surf_land 146')
+         enddo
+        enddo
+#endif
+#endif
 
 
@@ -172 +228 @@
        END DO
 
+#ifdef ISO
+      CALL abort_gcm('surf_land_mod 220','isos pas prevus dans orchidee',1)
+#endif
        ! temporary for keeping same results using lwdown_m instead of lwdown
        CALL surf_land_orchidee(itime, dtime, date0, knon, &
@@ -183 +242 @@
             tsol_rad, tsurf_new, alb1_new, alb2_new, &
             emis_new, z0m, z0h, qsurf, &
-            veget, lai, height)       
+            veget, lai, height &
+!#ifdef ISO
+!            , xtprecip_rain, xtprecip_snow, xtspechum, xtevap &
+!#endif
+            )                  
+
+#ifdef ISO
+#ifdef ISOVERIF
+     write(*,*) 'surf_land 193: apres surf_land_orchidee'   
+     do i=1,knon
+        if (iso_eau.gt.0) then
+             call iso_verif_egalite_choix(xtevap(iso_eau,i),evap(i), &
+    &            'surf_land 197',errmax,errmaxrel)
+        endif !if (iso_eau.gt.0) then      
+      enddo !do i=1,knon  
+#endif
+#endif
 !  
 !* Add contribution of relief to surface roughness
@@ -196 +271 @@
 !
 !****************************************************************************************
+#ifdef ISO
+#ifdef ISOVERIF
+!       write(*,*) 'surf_land 247'
+        call iso_verif_egalite_vect1D( &
+     &           xtsnow,snow,'surf_land_mod 207',niso,klon)
+#endif
+#endif
+
+#ifdef ISO
+        if (nudge_qsol.eq.1) then
+          call surf_land_nudge_qsol(knon,rlat,rlon,qsol,xtsol,knindex)
+        endif
+        !write(*,*) 'surf_land 258'
+#endif
        CALL surf_land_bucket(itime, jour, knon, knindex, debut, dtime,&
             tsurf, p1lay, cdragh, precip_rain, precip_totsnow, temp_air, &
@@ -202 +291 @@
             snow, qsol, agesno, tsoil, &
             qsurf, z0m, alb1_new, alb2_new, evap, &
-            fluxsens, fluxlat, tsurf_new, dflux_s, dflux_l)
+            fluxsens, fluxlat, tsurf_new, dflux_s, dflux_l &
+#ifdef ISO
+            ,xtprecip_rain, xtprecip_snow,xtspechum, &
+            xtsnow, xtsol,xtevap,h1, &
+     &      runoff_diag, xtrunoff_diag,Rland_ice &
+#endif           
+     &       )
         z0h(1:knon)=z0m(1:knon) ! En attendant mieux
 
@@ -224 +319 @@
          p1lay, temp_air, &
          flux_u1, flux_v1)
+
+#ifdef ISO
+#ifdef ISOVERIF
+!     write(*,*) 'surf_land 237: sortie'   
+      DO i=1,knon
+        IF (iso_eau >= 0) THEN
+             call iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
+    &            'surf_land 241',errmax,errmaxrel)
+        ENDIF !if (iso_eau.gt.0) then      
+      ENDDO !do i=1,knon  
+#endif
+#endif
 
 !albedo SB >>>
@@ -248 +355 @@
     
   END SUBROUTINE surf_land
+
+
+#ifdef ISO
+  SUBROUTINE surf_land_nudge_qsol(knon,rlat,rlon,qsol,xtsol,knindex)
+
+    USE dimphy   
+    USE infotrac_phy, ONLY: niso
+    USE isotopes_mod, ONLY: region_nudge_qsol    
+    INTEGER, INTENT(IN)                       :: knon         
+    REAL, DIMENSION(klon), INTENT(IN)         :: rlon, rlat
+    REAL, DIMENSION(klon), INTENT(INOUT)      :: qsol
+    INTEGER, DIMENSION(klon), INTENT(IN)      :: knindex    
+    REAL, DIMENSION(niso,klon), INTENT(INOUT) :: xtsol 
+    REAL :: lat_min_nudge_qsol,lat_max_nudge_qsol
+    REAL :: lon_min_nudge_qsol,lon_max_nudge_qsol
+    INTEGER :: i,ixt
+    REAL :: qsol_new
+
+    IF (region_nudge_qsol == 1) THEN
+        ! Aamzonie du Sud
+        lat_min_nudge_qsol=-15.0
+        lat_max_nudge_qsol=-5.0
+        lon_min_nudge_qsol=-70.0
+        lon_max_nudge_qsol=-50.0
+    ELSE IF (region_nudge_qsol == 2) THEN
+        ! Aamzonie du Nord
+        lat_min_nudge_qsol=-5.0
+        lat_max_nudge_qsol=5.0
+        lon_min_nudge_qsol=-70.0
+        lon_max_nudge_qsol=-50.0
+    ELSE
+        WRITE(*,*) 'surf_land 298: cas pas prevu'
+        WRITE(*,*) 'region_nudge_qsol=',region_nudge_qsol
+        stop
+    ENDIF
+
+!    write(*,*) 'surf_land 314: knon=',knon
+!    write(*,*) 'rlat=',rlat
+!    write(*,*) 'rlon=',rlon
+!    write(*,*) 'region_nudge_qsol=',region_nudge_qsol
+
+    DO i=1,knon
+      IF ((rlat(knindex(i)) >= lat_min_nudge_qsol).and. &
+  &       (rlat(knindex(i)) <= lat_max_nudge_qsol).and. &
+  &       (rlon(knindex(i)) >= lon_min_nudge_qsol).and. &
+  &       (rlon(knindex(i)) <= lon_max_nudge_qsol)) THEN
+!        write(*,*) 'surf_land 324: bon domaine: rlat,rlon,qsol=', &
+!  &             rlat(knindex(i)),rlon(knindex(i)),qsol(knindex(i))
+        qsol_new=qsol(i)
+        IF (region_nudge_qsol == 1) THEN   
+           qsol_new=max(qsol(i),50.0)   
+        ELSE IF (region_nudge_qsol == 2) THEN      
+           qsol_new=max(qsol(i),120.0)
+        ELSE !if (region_nudge_qsol.eq.1) then
+           WRITE(*,*) 'surf_land 317: cas pas prevu'
+           WRITE(*,*) 'region_nudge_qsol=',region_nudge_qsol
+           STOP
+        ENDIF !if (region_nudge_qsol.eq.1) then
+        IF (qsol(i) > 0.0) THEN
+           DO ixt=1,niso
+              xtsol(ixt,i)=xtsol(ixt,i)*qsol_new/qsol(i)
+           ENDDO
+        ELSE !IF (qsol(i) > 0.0) THEN
+           DO ixt=1,niso
+             xtsol(ixt,i)=0.0
+           ENDDO
+        ENDIF !IF (qsol(i) > 0.0) THEN
+        qsol(i)=qsol_new
+        WRITE(*,*) 'surf_land 346: qsol_new=',qsol(i)      
+     ENDIF ! if ((rlat(i).ge.lat_min_nudge_qsol).and.
+  ENDDO !DO i=1,knon
+
+  END SUBROUTINE surf_land_nudge_qsol
+#endif
+
 !
 !****************************************************************************************

LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90

@@ -23 +23 @@
        snowhgt, qsnow, to_ice, sissnow, &
        alb3, runoff, &
-       flux_u1, flux_v1)
+       flux_u1, flux_v1 &
+#ifdef ISO
+         &      ,xtprecip_rain, xtprecip_snow,xtspechum,Rland_ice &
+         &      ,xtsnow,xtsol,xtevap &
+#endif               
+           &    )
 
     USE dimphy
@@ -33 +38 @@
     USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic, zxqsaltlic
     USE phys_output_var_mod, ONLY : snow_o,zfra_o
+#ifdef ISO   
+    USE fonte_neige_mod,  ONLY : xtrun_off_lic
+    USE infotrac_phy,     ONLY : ntiso,niso
+    USE isotopes_routines_mod, ONLY: calcul_iso_surf_lic_vectall
+#ifdef ISOVERIF
+    USE isotopes_mod, ONLY: iso_eau,ridicule
+    USE isotopes_verif_mod
+#endif
+#endif
+ 
 !FC
     USE ioipsl_getin_p_mod, ONLY : getin_p
@@ -68 +83 @@
     REAL, DIMENSION(klon), INTENT(IN)             :: rugoro
     REAL, DIMENSION(klon,nbsrf), INTENT(IN)       :: pctsrf
+#ifdef ISO
+    REAL, DIMENSION(ntiso,klon), INTENT(IN)       :: xtprecip_rain, xtprecip_snow 
+    REAL, DIMENSION(ntiso,klon), INTENT(IN)       :: xtspechum 
+#endif
+
 
     LOGICAL,  INTENT(IN)                          :: debut   !true if first step
@@ -85 +105 @@
     REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
     REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
+#ifdef ISO
+    REAL, DIMENSION(niso,klon), INTENT(INOUT)     :: xtsnow, xtsol 
+    REAL, DIMENSION(niso,klon), INTENT(INOUT)     :: Rland_ice
+#endif
+
 
 ! Output variables
@@ -108 +133 @@
     REAL, DIMENSION(klon), INTENT(OUT)           :: sissnow
     REAL, DIMENSION(klon), INTENT(OUT)           :: runoff  !Land ice runoff
+#ifdef ISO
+    REAL, DIMENSION(ntiso,klon), INTENT(OUT)     :: xtevap     
+!    real, DIMENSION(niso,klon) :: xtrun_off_lic_0_diag ! est une variable globale de
+!    fonte_neige
+#endif
  
 
@@ -120 +150 @@
     REAL, DIMENSION(klon)    :: fqfonte,ffonte
     REAL, DIMENSION(klon)    :: run_off_lic_frac
+#ifdef ISO       
+    REAL, PARAMETER          :: t_coup = 273.15
+    REAL, DIMENSION(klon)    :: fqfonte_diag
+    REAL, DIMENSION(klon)    :: fq_fonte_diag
+    REAL, DIMENSION(klon)    ::  snow_evap_diag 
+    REAL, DIMENSION(klon)    ::  fqcalving_diag 
+    REAL max_eau_sol_diag  
+    REAL, DIMENSION(klon)    ::  runoff_diag 
+    REAL, DIMENSION(klon)    ::    run_off_lic_diag 
+    REAL                     ::  coeff_rel_diag
+    INTEGER                  :: ixt
+    REAL, DIMENSION(niso,klon) :: xtsnow_prec,xtsol_prec
+    REAL, DIMENSION(klon) :: snow_prec,qsol_prec
+!    real, DIMENSION(klon) :: run_off_lic_0_diag
+#endif
+
+
     REAL, DIMENSION(klon)    :: emis_new                  !Emissivity
     REAL, DIMENSION(klon)    :: swdown,lwdown
@@ -163 +210 @@
 !FC firtscall initializations
 !******************************************************************************************
+#ifdef ISO
+#ifdef ISOVERIF
+!     write(*,*) 'surf_land_ice 1499'   
+  DO i=1,knon
+    IF (iso_eau > 0) THEN
+      CALL iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
+    &                              'surf_land_ice 126',errmax,errmaxrel)
+    ENDIF !IF (iso_eau > 0) THEN      
+  ENDDO !DO i=1,knon  
+#endif
+#endif
+
   IF (firstcall) THEN
   alb_vis_sno_lic=0.77
@@ -202 +261 @@
 !****************************************************************************************
 #ifdef CPP_INLANDSIS
+
+#ifdef ISO
+        CALL abort_gcm('surf_landice 235','isotopes pas dans INLANDSIS',1)
+#endif
 
         debut_is=debut
@@ -322 +385 @@
          1.,AcoefH, AcoefQ, BcoefH, BcoefQ, &
          tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+
+#ifdef ISO
+#ifdef ISOVERIF
+     !write(*,*) 'surf_land_ice 1499'   
+     DO i=1,knon
+       IF (iso_eau > 0) THEN
+         IF (snow(i) > ridicule) THEN
+           CALL iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
+    &                                   'surf_land_ice 1151',errmax,errmaxrel)
+         ENDIF !IF ((snow(i) > ridicule)) THEN
+       ENDIF !IF (iso_eau > 0) THEN
+     ENDDO !DO i=1,knon  
+#endif
+
+    DO i=1,knon
+      snow_prec(i)=snow(i)
+      DO ixt=1,niso
+        xtsnow_prec(ixt,i)=xtsnow(ixt,i)
+      ENDDO !DO ixt=1,niso
+      ! initialisation:
+      fq_fonte_diag(i)=0.0
+      fqfonte_diag(i)=0.0
+      snow_evap_diag(i)=0.0 
+    ENDDO !DO i=1,knon 
+#endif         
 
     CALL calcul_flux_wind(knon, dtime, &
@@ -523 +611 @@
  
     CALL fonte_neige(knon, is_lic, knindex, dtime, &
-         tsurf, precip_rain, precip_totsnow,  &
-         snow, qsol, tsurf_new, evap_totsnow)
+         tsurf, precip_rain, precip_totsnow, &
+         snow, qsol, tsurf_new, evap_totsnow &
+#ifdef ISO    
+     &  ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag     &
+     &  ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag &
+#endif
+     &   )
+
+
+#ifdef ISO
+#ifdef ISOVERIF
+    DO i=1,knon  
+      IF (iso_eau > 0) THEN  
+        CALL iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, &
+     &                               'surf_landice_mod 217',errmax,errmaxrel)
+      ENDIF !IF (iso_eau > 0) THEN
+    ENDDO !DO i=1,knon
+#endif
+
+    CALL calcul_iso_surf_lic_vectall(klon,knon, &
+     &    evap,snow_evap_diag,Tsurf_new,snow, &
+     &    fq_fonte_diag,fqfonte_diag,dtime,t_coup, &
+     &    precip_snow,xtprecip_snow,precip_rain,xtprecip_rain, snow_prec,xtsnow_prec, &
+     &    xtspechum,spechum,ps,Rland_ice, &
+     &    xtevap,xtsnow,fqcalving_diag, &
+     &    knindex,is_lic,run_off_lic_diag,coeff_rel_diag &
+     &   ) 
+
+!        call fonte_neige_export_xtrun_off_lic_0(knon,xtrun_off_lic_0_diag)
+
+#endif
    
-    
     WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.                                         
     zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))  

LMDZ6/trunk/libf/phylmd/surf_ocean_mod.F90

@@ -21 +21 @@
        tsurf_new, dflux_s, dflux_l, lmt_bils, &
        flux_u1, flux_v1, delta_sst, delta_sal, ds_ns, dt_ns, dter, dser, &
-       dt_ds, tkt, tks, taur, sss)
+       dt_ds, tkt, tks, taur, sss &
+#ifdef ISO
+        &       ,xtprecip_rain, xtprecip_snow,xtspechum,Roce, &
+        &       xtsnow,xtevap,h1 &
+#endif               
+        &       )
 
     use albedo, only: alboc, alboc_cd
@@ -31 +36 @@
     USE ocean_cpl_mod, ONLY    : ocean_cpl_noice
     USE indice_sol_mod, ONLY : nbsrf, is_oce
+#ifdef ISO
+    USE infotrac_phy, ONLY : ntraciso=>ntiso,niso
+#ifdef ISOVERIF
+    USE isotopes_mod, ONLY: iso_eau,ridicule
+    USE isotopes_verif_mod
+#endif
+#endif
     USE limit_read_mod
-    use config_ocean_skin_m, only: activate_ocean_skin
+    USE config_ocean_skin_m, ONLY: activate_ocean_skin
     !
     ! This subroutine will make a call to ocean_XXX_noice according to the ocean mode (force, 
@@ -68 +80 @@
     REAL, DIMENSION(klon), INTENT(IN)        :: rugoro
     REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow 
+    REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: xtspechum
+#endif
 
     ! In/Output variables
@@ -75 +91 @@
     REAL, DIMENSION(klon), INTENT(INOUT)     :: agesno
     REAL, DIMENSION(klon), INTENT(inOUT)     :: z0h
+#ifdef ISO
+    REAL, DIMENSION(niso,klon), INTENT(IN)   :: xtsnow 
+    REAL, DIMENSION(niso,klon), INTENT(INOUT):: Roce  
+#endif
 
     REAL, intent(inout):: delta_sst(:) ! (knon)
@@ -136 +156 @@
     ! size klon because of the coupling machinery.)
 
+#ifdef ISO
+    REAL, DIMENSION(ntraciso,klon), INTENT(out) :: xtevap ! isotopes in surface evaporation flux 
+    REAL, DIMENSION(klon), INTENT(out)          :: h1 ! just a diagnostic, not useful for the simulation   
+#endif
+
     ! Local variables
     !*************************************************************************
@@ -146 +171 @@
     REAL, DIMENSION(klon) :: precip_totsnow
     CHARACTER(len=20),PARAMETER :: modname="surf_ocean"
-    real rhoa(knon) ! density of moist air  (kg / m3)
+    REAL rhoa(knon) ! density of moist air  (kg / m3)
     REAL sens_prec_liq(knon)
 
     REAL t_int(knon) ! ocean-air interface temperature, in K
-    real s_int(knon) ! ocean-air interface salinity, in ppt
+    REAL s_int(knon) ! ocean-air interface salinity, in ppt
 
     !**************************************************************************
 
+#ifdef ISO
+#ifdef ISOVERIF
+    DO i = 1, knon
+      IF (iso_eau > 0) THEN         
+        CALL iso_verif_egalite_choix(xtspechum(iso_eau,i), &
+     &          spechum(i),'surf_ocean_mod 117', &
+     &          errmax,errmaxrel)          
+        CALL iso_verif_egalite_choix(xtsnow(iso_eau,i), &
+     &          snow(i),'surf_ocean_mod 127', &
+     &          errmax,errmaxrel) 
+      ENDIF !IF (iso_eau > 0) then
+    ENDDO !DO i=1,klon
+#endif      
+#endif
 
     !******************************************************************************
@@ -230 +269 @@
             radsol, snow, agesno, &
             qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-            tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa)
+            tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa &
+#ifdef ISO
+            ,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, &
+            xtsnow,xtevap,h1 &  
+#endif            
+            )
     END SELECT
 

LMDZ6/trunk/libf/phylmd/surf_seaice_mod.F90

@@ -21 +21 @@
        z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, &  
        tsurf_new, dflux_s, dflux_l, &
-       flux_u1, flux_v1)
+       flux_u1, flux_v1 &
+#ifdef ISO
+         &      ,xtprecip_rain, xtprecip_snow,xtspechum,Roce, &
+         &      xtsnow,xtsol,xtevap,Rland_ice &
+#endif               
+         &      )
 
   USE dimphy
@@ -29 +34 @@
   USE ocean_slab_mod, ONLY   : ocean_slab_ice
   USE indice_sol_mod
+#ifdef ISO
+  USE infotrac_phy, ONLY : ntiso,niso
+#endif
 
 !
@@ -62 +70 @@
     REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1, gustiness
     REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf
+#ifdef ISO
+    REAL, DIMENSION(ntiso,klon), INTENT(IN)  :: xtprecip_rain, xtprecip_snow 
+    REAL, DIMENSION(klon),       INTENT(IN)  :: xtspechum
+    REAL, DIMENSION(niso,klon),  INTENT(IN)  :: Roce
+    REAL, DIMENSION(niso,klon),  INTENT(IN)  :: Rland_ice
+#endif
 
 ! In/Output arguments
@@ -68 +82 @@
     REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
     REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
+#ifdef ISO
+    REAL, DIMENSION(niso,klon), INTENT(INOUT)     :: xtsnow  
+    REAL, DIMENSION(niso,klon), INTENT(IN)        :: xtsol 
+#endif
 
 ! Output arguments
@@ -82 +100 @@
     REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
     REAL, DIMENSION(klon), INTENT(OUT)       :: flux_u1, flux_v1
+#ifdef ISO
+    REAL, DIMENSION(ntiso,klon), INTENT(OUT) :: xtevap
+#endif
 
 ! Local arguments
 !****************************************************************************************
     REAL, DIMENSION(klon)  :: radsol
+#ifdef ISO
+#ifdef ISOVERIF
+    INTEGER :: j
+#endif
+#endif
 
 !albedo SB >>>
@@ -145 +171 @@
             radsol, snow, qsol, agesno, tsoil, &
             qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-            tsurf_new, dflux_s, dflux_l, rhoa)
+            tsurf_new, dflux_s, dflux_l, rhoa &
+#ifdef ISO
+            ,xtprecip_rain, xtprecip_snow, xtspechum,Roce, &
+            xtsnow, xtsol,xtevap,Rland_ice &  
+#endif            
+            )
 
     END IF