Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/aeropt_2bands.F90
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/aeropt_2bands.F90	(revision 1346)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/aeropt_2bands.F90	(revision 1347)
@@ -10,4 +10,5 @@
   USE dimphy
   USE aero_mod
+  USE phys_local_var_mod, only: absvisaer
 
   !    Yves Balkanski le 12 avril 2006
@@ -1119,4 +1120,10 @@
     ENDDO
   ENDDO
+   
+
+  inu=1
+  DO i=1, KLON
+     absvisaer(i)=SUM((1-piz_allaer(i,:,:,inu))*tau_allaer(i,:,:,inu))
+  END DO	
 
   DEALLOCATE(aerosol_name) 
Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/aeropt_5wv.F90
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/aeropt_5wv.F90	(revision 1346)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/aeropt_5wv.F90	(revision 1347)
@@ -11,4 +11,5 @@
   USE DIMPHY
   USE aero_mod
+  USE phys_local_var_mod, only: od550aer,od865aer,ec550aer,od550lt1aer
 
   !
@@ -141,5 +142,5 @@
   REAL :: tau3d(KLON,KLEV), piz3d(KLON,KLEV), cg3d(KLON,KLEV)
   REAL :: abs3d(KLON,KLEV)     ! epaisseur optique d'absorption
-
+  REAL :: dh(KLON,KLEV)
   
   REAL :: alpha_aers_5wv(nbre_RH,las,naero_soluble)   ! ext. coeff. Soluble comp. units *** m2/g 
@@ -621,5 +622,8 @@
   DO k=1, klev
     DO i=1, klon
+!      IF (t_seri(i,k).EQ.0) stop 'stop aeropt_5wv T '
+!      IF (pplay(i,k).EQ.0) stop  'stop aeropt_5wv p '
       zrho=pplay(i,k)/t_seri(i,k)/RD                  ! kg/m3
+      dh(i,k)=pdel(i,k)/(gravit*zrho)
 !CDIR UNROLL=naero_spc
       mass_temp(i,k,:) = m_allaer(i,k,:) / zrho / 1.e+9
@@ -846,4 +850,17 @@
 !  ENDDO
 
+   DO i=1, klon
+      od550aer(i)=SUM(tausum(i,2,:))
+      od865aer(i)=SUM(tausum(i,5,:))
+      DO k=1, KLEV
+         ec550aer(i,k)=SUM(tau(i,k,2,:))/dh(i,k)
+      END DO	
+   END DO
+   od550lt1aer(:)=tausum(:,2,id_ASSO4M)+tausum(:,2,id_ASBCM)+tausum(:,2,id_AIBCM)+ &
+	tausum(:,2,id_ASPOMM)+tausum(:,2,id_AIPOMM)+tausum(:,2,id_ASSSM)+ &
+	0.03*tausum(:,2,id_CSSSM)+0.4*tausum(:,2,id_CIDUSTM)
+
+
+
   DEALLOCATE(aerosol_name) 
   
Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/newmicro.F
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/newmicro.F	(revision 1346)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/newmicro.F	(revision 1347)
@@ -11,4 +11,8 @@
 
       USE dimphy
+      USE phys_local_var_mod, only: scdnc,cldncl,reffclwtop,lcc,
+     .                              reffclws,reffclwc,cldnvi,lcc3d,
+     .                              lcc3dcon,lcc3dstra
+      USE phys_state_var_mod, only: rnebcon,clwcon
       IMPLICIT none
 c======================================================================
@@ -46,4 +50,19 @@
 #include "radepsi.h"
 #include "radopt.h"
+c choix de l'hypothese de recouvrememnt nuageuse
+      LOGICAL RANDOM,MAXIMUM_RANDOM,MAXIMUM,FIRST
+      parameter (RANDOM=.FALSE., MAXIMUM_RANDOM=.TRUE., MAXIMUM=.FALSE.)
+c Hypoyhese de recouvrement : MAXIMUM_RANDOM
+      INTEGER flag_max
+      REAL phase3d(klon, klev),dh(klon, klev),pdel(klon, klev),
+     .     zrho(klon, klev)
+      REAL tcc(klon), ftmp(klon), lcc_integrat(klon), height(klon)
+      REAL thres_tau,thres_neb
+      PARAMETER (thres_tau=0.3, thres_neb=0.001)
+      REAL t_tmp
+      REAL gravit
+      PARAMETER (gravit=9.80616)  !m/s2
+      REAL pqlwpcon(klon, klev), pqlwpstra(klon, klev) 
+c
       REAL paprs(klon,klev+1), pplay(klon,klev)
       REAL t(klon,klev)
@@ -131,4 +150,8 @@
       xflwc = 0.d0
       xfiwc = 0.d0
+
+! Initialisation 
+      reliq=0.
+      reice=0.
 
       DO k = 1, klev
@@ -471,5 +494,183 @@
          pcl(i)=1.-pcl(i)
       ENDDO
-      
+
+c ========================================================
+! DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL
+c ========================================================
+!! change by Nicolas Yan (LSCE) 
+!! Cloud Droplet Number Concentration (CDNC) : 3D variable
+!! Fractionnal cover by liquid water cloud (LCC3D) : 3D variable
+!! Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable
+!! Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable
+!! Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable
+      IF (ok_newmicro) THEN
+         IF (ok_aie) THEN
+            DO k = 1, klev
+               DO i = 1, klon
+                  phase3d(i,k)=1-zfice2(i,k)
+                  IF (pclc(i,k) .LE. seuil_neb) THEN
+                     lcc3d(i,k)=seuil_neb*phase3d(i,k)
+                  ELSE
+                     lcc3d(i,k)=pclc(i,k)*phase3d(i,k)
+                  ENDIF
+                  scdnc(i,k)=lcc3d(i,k)*cdnc(i,k) ! m-3
+               ENDDO
+            ENDDO
+
+            DO i=1,klon
+               lcc(i)=0.
+               reffclwtop(i)=0.
+               cldncl(i)=0.
+               IF(RANDOM .OR. MAXIMUM_RANDOM) tcc(i) = 1.
+               IF(MAXIMUM) tcc(i) = 0.
+            ENDDO
+     
+            FIRST=.TRUE.
+
+            DO i=1,klon
+               DO k=klev-1,1,-1 !From TOA down
+
+
+            ! Test, if the cloud optical depth exceeds the necessary
+            ! threshold:
+
+             IF (pcltau(i,k).GT.thres_tau .AND. pclc(i,k).GT.thres_neb)
+     .                                                             THEN
+               ! To calculate the right Temperature at cloud top,
+               ! interpolate it between layers:
+                  t_tmp = t(i,k) +
+     .              (paprs(i,k+1)-pplay(i,k))/(pplay(i,k+1)-pplay(i,k))
+     .              * ( t(i,k+1) - t(i,k) )
+
+                  IF(MAXIMUM) THEN
+                    IF(FIRST) THEN
+                       write(*,*)'Hypothese de recouvrement: MAXIMUM'
+                       FIRST=.FALSE.
+                    ENDIF
+                    flag_max= -1.
+                    ftmp(i) = MAX(tcc(i),pclc(i,k))
+                  ENDIF
+
+                  IF(RANDOM) THEN
+                    IF(FIRST) THEN
+                       write(*,*)'Hypothese de recouvrement: RANDOM'
+                       FIRST=.FALSE.
+                    ENDIF
+                    flag_max= 1.
+                    ftmp(i) = tcc(i) * (1-pclc(i,k))
+                  ENDIF
+
+                  IF(MAXIMUM_RANDOM) THEN
+                    IF(FIRST) THEN
+                       write(*,*)'Hypothese de recouvrement: MAXIMUM_
+     .                         RANDOM'
+                       FIRST=.FALSE.
+                    ENDIF
+                    flag_max= 1.
+                    ftmp(i) = tcc(i) *
+     .              (1. - MAX(pclc(i,k),pclc(i,k+1))) /
+     .              (1. - MIN(pclc(i,k+1),1.-thres_neb))
+                  ENDIF
+c Effective radius of cloud droplet at top of cloud (m)
+                  reffclwtop(i) = reffclwtop(i) + rad_chaud_tab(i,k) * 
+     .           1.0E-06 * phase3d(i,k) * ( tcc(i) - ftmp(i))*flag_max
+c CDNC at top of cloud (m-3)
+                  cldncl(i) = cldncl(i) + cdnc(i,k) * phase3d(i,k) * 
+     .                 (tcc(i) - ftmp(i))*flag_max
+c Liquid Cloud Content at top of cloud
+                  lcc(i) = lcc(i) + phase3d(i,k) * (tcc(i)-ftmp(i))*
+     .                    flag_max
+c Total Cloud Content at top of cloud
+                  tcc(i)=ftmp(i)
+              
+          ENDIF ! is there a visible, not-too-small cloud?  
+          ENDDO ! loop over k
+
+          IF(RANDOM .OR. MAXIMUM_RANDOM) tcc(i)=1.-tcc(i)
+         ENDDO ! loop over i
+
+!! Convective and Stratiform Cloud Droplet Effective Radius (REFFCLWC  REFFCLWS)
+            DO i = 1, klon
+               DO k = 1, klev
+                  pqlwpcon(i,k)=rnebcon(i,k)*clwcon(i,k) ! fraction eau liquide convective
+                  pqlwpstra(i,k)=pclc(i,k)*phase3d(i,k)-pqlwpcon(i,k) ! fraction eau liquide stratiforme
+                  IF (pqlwpstra(i,k) .LE. 0.0) pqlwpstra(i,k)=0.0
+! Convective Cloud Droplet Effective Radius (REFFCLWC) : variable 3D
+                  reffclwc(i,k)=1.1
+     &                 *((pqlwpcon(i,k)*pplay(i,k)/(RD * T(i,k)))
+     &                 /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.)
+                  reffclwc(i,k) = MAX(reffclwc(i,k) * 1e6, 5.)
+
+! Stratiform Cloud Droplet Effective Radius (REFFCLWS) : variable 3D
+                  IF ((pclc(i,k)-rnebcon(i,k)) .LE. seuil_neb) THEN ! tout sous la forme convective
+                     reffclws(i,k)=0.0
+                     lcc3dstra(i,k)= 0.0
+                  ELSE
+                     reffclws(i,k) = (pclc(i,k)*phase3d(i,k)*
+     &                               rad_chaud_tab(i,k)-
+     &                            pqlwpcon(i,k)*reffclwc(i,k))
+                     IF(reffclws(i,k) .LE. 0.0) reffclws(i,k)=0.0
+                     lcc3dstra(i,k)=pqlwpstra(i,k)
+                 ENDIF
+!Convertion from um to m
+                  IF(rnebcon(i,k). LE. seuil_neb) THEN
+                    reffclwc(i,k) = reffclwc(i,k)*seuil_neb*clwcon(i,k)
+     &                              *1.0E-06
+                    lcc3dcon(i,k)= seuil_neb*clwcon(i,k)
+                  ELSE
+                    reffclwc(i,k) = reffclwc(i,k)*pqlwpcon(i,k)
+     &                              *1.0E-06
+                    lcc3dcon(i,k) = pqlwpcon(i,k)
+                  ENDIF
+
+                  reffclws(i,k) = reffclws(i,k)*1.0E-06
+
+               ENDDO !klev
+            ENDDO !klon 
+
+!! Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D
+            DO k = 1, klev
+               DO i = 1, klon
+                   pdel(i,k) = paprs(i,k)-paprs(i,k+1)
+                   zrho(i,k)=pplay(i,k)/t(i,k)/RD                  ! kg/m3
+                   dh(i,k)=pdel(i,k)/(gravit*zrho(i,k)) ! hauteur de chaque boite (m)
+               ENDDO
+            ENDDO
+c
+            DO i = 1, klon
+               cldnvi(i)=0.
+               lcc_integrat(i)=0.
+               height(i)=0.
+               DO k = 1, klev
+                  cldnvi(i)=cldnvi(i)+cdnc(i,k)*lcc3d(i,k)*dh(i,k)
+                  lcc_integrat(i)=lcc_integrat(i)+lcc3d(i,k)*dh(i,k)
+                  height(i)=height(i)+dh(i,k)
+               ENDDO ! klev
+               lcc_integrat(i)=lcc_integrat(i)/height(i)
+               IF (lcc_integrat(i) .LE. 1.0E-03) THEN
+                  cldnvi(i)=cldnvi(i)*lcc(i)/seuil_neb
+               ELSE
+                  cldnvi(i)=cldnvi(i)*lcc(i)/lcc_integrat(i)
+               ENDIF
+            ENDDO ! klon
+            
+            DO i = 1, klon
+               DO k = 1, klev
+                IF (scdnc(i,k) .LE. 0.0) scdnc(i,k)=0.0
+                IF (reffclws(i,k) .LE. 0.0) reffclws(i,k)=0.0
+                IF (reffclwc(i,k) .LE. 0.0) reffclwc(i,k)=0.0
+                IF (lcc3d(i,k) .LE. 0.0) lcc3d(i,k)=0.0
+                IF (lcc3dcon(i,k) .LE. 0.0) lcc3dcon(i,k)=0.0
+                IF (lcc3dstra(i,k) .LE. 0.0) lcc3dstra(i,k)=0.0
+               ENDDO
+               IF (reffclwtop(i) .LE. 0.0) reffclwtop(i)=0.0
+               IF (cldncl(i) .LE. 0.0) cldncl(i)=0.0
+               IF (cldnvi(i) .LE. 0.0) cldnvi(i)=0.0
+               IF (lcc(i) .LE. 0.0) lcc(i)=0.0
+            ENDDO
+
+         ENDIF !ok_aie
+      ENDIF !ok newmicro
+c
 C
       RETURN
Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phys_local_var_mod.F90	(revision 1346)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phys_local_var_mod.F90	(revision 1347)
@@ -88,4 +88,78 @@
       REAL, SAVE, ALLOCATABLE :: tau3d_aero(:,:,:,:) 
       !$OMP THREADPRIVATE(tau3d_aero) 
+      REAL, SAVE, ALLOCATABLE :: scdnc(:,:)
+      !$OMP THREADPRIVATE(scdnc)
+      REAL, SAVE, ALLOCATABLE :: cldncl(:)
+      !$OMP THREADPRIVATE(cldncl)
+      REAL, SAVE, ALLOCATABLE :: reffclwtop(:)
+      !$OMP THREADPRIVATE(reffclwtop)
+      REAL, SAVE, ALLOCATABLE :: lcc(:)
+      !$OMP THREADPRIVATE(lcc)
+      REAL, SAVE, ALLOCATABLE :: reffclws(:,:)
+      !$OMP THREADPRIVATE(reffclws)
+      REAL, SAVE, ALLOCATABLE :: reffclwc(:,:)
+      !$OMP THREADPRIVATE(reffclwc)
+      REAL, SAVE, ALLOCATABLE :: cldnvi(:) 
+      !$OMP THREADPRIVATE(cldnvi)
+      REAL, SAVE, ALLOCATABLE :: lcc3d(:,:)
+      !$OMP THREADPRIVATE(lcc3d)
+      REAL, SAVE, ALLOCATABLE :: lcc3dcon(:,:)
+      !$OMP THREADPRIVATE(lcc3dcon)
+      REAL, SAVE, ALLOCATABLE :: lcc3dstra(:,:)
+      !$OMP THREADPRIVATE(lcc3dstra)
+      REAL, SAVE, ALLOCATABLE :: od550aer(:) 
+      !$OMP THREADPRIVATE(od550aer) 
+      REAL, SAVE, ALLOCATABLE :: absvisaer(:) 
+      !$OMP THREADPRIVATE(absvisaer) 
+      REAL, SAVE, ALLOCATABLE :: od865aer(:) 
+      !$OMP THREADPRIVATE(od865aer) 
+      REAL, SAVE, ALLOCATABLE :: ec550aer(:,:) 
+      !$OMP THREADPRIVATE(ec550aer) 
+      REAL, SAVE, ALLOCATABLE :: od550lt1aer(:) 
+      !$OMP THREADPRIVATE(od550lt1aer) 
+      REAL, SAVE, ALLOCATABLE :: sconcso4(:) 
+      !$OMP THREADPRIVATE(sconcso4) 
+      REAL, SAVE, ALLOCATABLE :: sconcoa(:) 
+      !$OMP THREADPRIVATE(sconcoa) 
+      REAL, SAVE, ALLOCATABLE :: sconcbc(:) 
+      !$OMP THREADPRIVATE(sconcbc) 
+      REAL, SAVE, ALLOCATABLE :: sconcss(:) 
+      !$OMP THREADPRIVATE(sconcss) 
+      REAL, SAVE, ALLOCATABLE :: sconcdust(:) 
+      !$OMP THREADPRIVATE(sconcdust) 
+      REAL, SAVE, ALLOCATABLE :: concso4(:,:) 
+      !$OMP THREADPRIVATE(concso4) 
+      REAL, SAVE, ALLOCATABLE :: concoa(:,:) 
+      !$OMP THREADPRIVATE(concoa) 
+      REAL, SAVE, ALLOCATABLE :: concbc(:,:) 
+      !$OMP THREADPRIVATE(concbc) 
+      REAL, SAVE, ALLOCATABLE :: concss(:,:) 
+      !$OMP THREADPRIVATE(concss) 
+      REAL, SAVE, ALLOCATABLE :: concdust(:,:) 
+      !$OMP THREADPRIVATE(concdust) 
+      REAL, SAVE, ALLOCATABLE :: loadso4(:) 
+      !$OMP THREADPRIVATE(loadso4) 
+      REAL, SAVE, ALLOCATABLE :: loadoa(:) 
+      !$OMP THREADPRIVATE(loadoa) 
+      REAL, SAVE, ALLOCATABLE :: loadbc(:) 
+      !$OMP THREADPRIVATE(loadbc) 
+      REAL, SAVE, ALLOCATABLE :: loadss(:) 
+      !$OMP THREADPRIVATE(loadss) 
+      REAL, SAVE, ALLOCATABLE :: loaddust(:) 
+      !$OMP THREADPRIVATE(loaddust) 
+      REAL, SAVE, ALLOCATABLE :: load_tmp1(:) 
+      !$OMP THREADPRIVATE(load_tmp1) 
+      REAL, SAVE, ALLOCATABLE :: load_tmp2(:) 
+      !$OMP THREADPRIVATE(load_tmp2) 
+      REAL, SAVE, ALLOCATABLE :: load_tmp3(:) 
+      !$OMP THREADPRIVATE(load_tmp3) 
+      REAL, SAVE, ALLOCATABLE :: load_tmp4(:) 
+      !$OMP THREADPRIVATE(load_tmp4) 
+      REAL, SAVE, ALLOCATABLE :: load_tmp5(:) 
+      !$OMP THREADPRIVATE(load_tmp5) 
+      REAL, SAVE, ALLOCATABLE :: load_tmp6(:) 
+      !$OMP THREADPRIVATE(load_tmp6) 
+      REAL, SAVE, ALLOCATABLE :: load_tmp7(:) 
+      !$OMP THREADPRIVATE(load_tmp7) 
 
 CONTAINS
@@ -131,4 +205,41 @@
       allocate(tausum_aero(klon,nwave,naero_spc))
       allocate(tau3d_aero(klon,klev,nwave,naero_spc)) 
+      allocate(scdnc(klon, klev))
+      allocate(cldncl(klon))
+      allocate(reffclwtop(klon))
+      allocate(lcc(klon))
+      allocate(reffclws(klon, klev))
+      allocate(reffclwc(klon, klev))
+      allocate(cldnvi(klon))
+      allocate(lcc3d(klon, klev))
+      allocate(lcc3dcon(klon, klev))
+      allocate(lcc3dstra(klon, klev))
+      allocate(od550aer(klon))	 
+      allocate(od865aer(klon))	 
+      allocate(absvisaer(klon))	 
+      allocate(ec550aer(klon,klev))
+      allocate(od550lt1aer(klon))	 	 
+      allocate(sconcso4(klon))
+      allocate(sconcoa(klon))
+      allocate(sconcbc(klon))
+      allocate(sconcss(klon))
+      allocate(sconcdust(klon))
+      allocate(concso4(klon,klev))
+      allocate(concoa(klon,klev))
+      allocate(concbc(klon,klev))
+      allocate(concss(klon,klev))
+      allocate(concdust(klon,klev))
+      allocate(loadso4(klon))
+      allocate(loadoa(klon))
+      allocate(loadbc(klon))
+      allocate(loadss(klon))
+      allocate(loaddust(klon))
+      allocate(load_tmp1(klon))
+      allocate(load_tmp2(klon))
+      allocate(load_tmp3(klon))
+      allocate(load_tmp4(klon))
+      allocate(load_tmp5(klon))
+      allocate(load_tmp6(klon))
+      allocate(load_tmp7(klon))
 
 END SUBROUTINE phys_local_var_init
@@ -170,4 +281,41 @@
       deallocate(tausum_aero) 
       deallocate(tau3d_aero) 
+      deallocate(scdnc)
+      deallocate(cldncl)
+      deallocate(reffclwtop)
+      deallocate(lcc)
+      deallocate(reffclws)
+      deallocate(reffclwc)
+      deallocate(cldnvi)
+      deallocate(lcc3d)
+      deallocate(lcc3dcon)
+      deallocate(lcc3dstra)
+      deallocate(od550aer)	 
+      deallocate(od865aer)
+      deallocate(absvisaer)
+      deallocate(ec550aer)
+      deallocate(od550lt1aer)
+      deallocate(sconcso4) 
+      deallocate(sconcoa) 
+      deallocate(sconcbc) 
+      deallocate(sconcss) 
+      deallocate(sconcdust) 
+      deallocate(concso4) 
+      deallocate(concoa) 
+      deallocate(concbc) 
+      deallocate(concss) 
+      deallocate(concdust) 
+      deallocate(loadso4) 
+      deallocate(loadoa) 
+      deallocate(loadbc) 
+      deallocate(loadss) 
+      deallocate(loaddust) 
+      deallocate(load_tmp1)
+      deallocate(load_tmp2)
+      deallocate(load_tmp3)
+      deallocate(load_tmp4)
+      deallocate(load_tmp5)
+      deallocate(load_tmp6)
+      deallocate(load_tmp7)
       deallocate(d_u_hin,d_v_hin,d_t_hin)
 
Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phys_output_mod.F90
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phys_output_mod.F90	(revision 1346)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phys_output_mod.F90	(revision 1347)
@@ -54,6 +54,6 @@
 
 !!! 1D
-  type(ctrl_out),save :: o_phis         = ctrl_out((/ 1, 1, 10, 1, 1 /), 'phis') 
-  type(ctrl_out),save :: o_aire         = ctrl_out((/ 1, 1, 10,  1, 1 /),'aire')
+  type(ctrl_out),save :: o_phis         = ctrl_out((/ 1, 1, 10, 5, 1 /), 'phis') 
+  type(ctrl_out),save :: o_aire         = ctrl_out((/ 1, 1, 10,  10, 1 /),'aire')
   type(ctrl_out),save :: o_contfracATM  = ctrl_out((/ 10, 1,  1, 10, 10 /),'contfracATM')
   type(ctrl_out),save :: o_contfracOR   = ctrl_out((/ 10, 1,  1, 10, 10 /),'contfracOR')
@@ -61,8 +61,8 @@
   
 !!! 2D
-  type(ctrl_out),save :: o_flat         = ctrl_out((/ 10, 1, 10, 10, 1 /),'flat')
+  type(ctrl_out),save :: o_flat         = ctrl_out((/ 5, 1, 10, 5, 1 /),'flat')
   type(ctrl_out),save :: o_slp          = ctrl_out((/ 1, 1, 1, 10, 1 /),'slp')
-  type(ctrl_out),save :: o_tsol         = ctrl_out((/ 1, 1, 1, 1, 1 /),'tsol')
-  type(ctrl_out),save :: o_t2m          = ctrl_out((/ 1, 1, 1, 1, 1 /),'t2m')
+  type(ctrl_out),save :: o_tsol         = ctrl_out((/ 1, 1, 1, 5, 1 /),'tsol')
+  type(ctrl_out),save :: o_t2m          = ctrl_out((/ 1, 1, 1, 5, 1 /),'t2m')
   type(ctrl_out),save :: o_t2m_min      = ctrl_out((/ 1, 1, 10, 10, 10 /),'t2m_min')
   type(ctrl_out),save :: o_t2m_max      = ctrl_out((/ 1, 1, 10, 10, 10 /),'t2m_max')
@@ -75,8 +75,8 @@
   type(ctrl_out),save :: o_wind10max    = ctrl_out((/ 10, 1, 10, 10, 10 /),'wind10max')
   type(ctrl_out),save :: o_sicf         = ctrl_out((/ 1, 1, 10, 10, 10 /),'sicf')
-  type(ctrl_out),save :: o_q2m          = ctrl_out((/ 1, 1, 1, 1, 1 /),'q2m')
-  type(ctrl_out),save :: o_u10m         = ctrl_out((/ 1, 1, 1, 1, 1 /),'u10m')
-  type(ctrl_out),save :: o_v10m         = ctrl_out((/ 1, 1, 1, 1, 1 /),'v10m')
-  type(ctrl_out),save :: o_psol         = ctrl_out((/ 1, 1, 1, 1, 1 /),'psol')
+  type(ctrl_out),save :: o_q2m          = ctrl_out((/ 1, 1, 1, 5, 1 /),'q2m')
+  type(ctrl_out),save :: o_u10m         = ctrl_out((/ 1, 1, 1, 5, 1 /),'u10m')
+  type(ctrl_out),save :: o_v10m         = ctrl_out((/ 1, 1, 1, 5, 1 /),'v10m')
+  type(ctrl_out),save :: o_psol         = ctrl_out((/ 1, 1, 1, 5, 1 /),'psol')
   type(ctrl_out),save :: o_qsurf        = ctrl_out((/ 1, 10, 10, 10, 10 /),'qsurf')
 
@@ -91,16 +91,23 @@
                                               ctrl_out((/ 10, 4, 10, 10, 10 /),'v10m_sic') /)
 
-  type(ctrl_out),save :: o_qsol         = ctrl_out((/ 1, 10, 10, 1, 1 /),'qsol')
+  type(ctrl_out),save :: o_qsol         = ctrl_out((/ 1, 10, 10, 10, 10 /),'qsol')
 
   type(ctrl_out),save :: o_ndayrain     = ctrl_out((/ 1, 10, 10, 10, 10 /),'ndayrain')
-  type(ctrl_out),save :: o_precip       = ctrl_out((/ 1, 1, 1, 1, 1 /),'precip')
-  type(ctrl_out),save :: o_plul         = ctrl_out((/ 1, 1, 1, 1, 10 /),'plul')
-
-  type(ctrl_out),save :: o_pluc         = ctrl_out((/ 1, 1, 1, 1, 10 /),'pluc')
-  type(ctrl_out),save :: o_snow         = ctrl_out((/ 1, 1, 10, 1, 10 /),'snow') 
-  type(ctrl_out),save :: o_evap         = ctrl_out((/ 1, 1, 10, 1, 10 /),'evap')
+  type(ctrl_out),save :: o_precip       = ctrl_out((/ 1, 1, 1, 5, 10 /),'precip')
+  type(ctrl_out),save :: o_plul         = ctrl_out((/ 1, 1, 1, 10, 10 /),'plul')
+
+  type(ctrl_out),save :: o_pluc         = ctrl_out((/ 1, 1, 1, 5, 10 /),'pluc')
+  type(ctrl_out),save :: o_snow         = ctrl_out((/ 1, 1, 10, 5, 10 /),'snow') 
+  type(ctrl_out),save :: o_evap         = ctrl_out((/ 1, 1, 10, 10, 10 /),'evap')
+  type(ctrl_out),save,dimension(4) :: o_evap_srf     = (/ ctrl_out((/ 1, 1, 10, 10, 10 /),'evap_ter'), &
+                                           ctrl_out((/ 1, 1, 10, 10, 10 /),'evap_lic'), &
+                                           ctrl_out((/ 1, 1, 10, 10, 10 /),'evap_oce'), &
+                                           ctrl_out((/ 1, 1, 10, 10, 10 /),'evap_sic') /)
+  type(ctrl_out),save :: o_msnow       = ctrl_out((/ 1, 10, 10, 10, 10 /),'msnow')
+  type(ctrl_out),save :: o_fsnow       = ctrl_out((/ 1, 10, 10, 10, 10 /),'fsnow')
+
   type(ctrl_out),save :: o_tops         = ctrl_out((/ 1, 1, 10, 10, 10 /),'tops')
   type(ctrl_out),save :: o_tops0        = ctrl_out((/ 1, 5, 10, 10, 10 /),'tops0')
-  type(ctrl_out),save :: o_topl         = ctrl_out((/ 1, 1, 10, 1, 10 /),'topl')
+  type(ctrl_out),save :: o_topl         = ctrl_out((/ 1, 1, 10, 5, 10 /),'topl')
   type(ctrl_out),save :: o_topl0        = ctrl_out((/ 1, 5, 10, 10, 10 /),'topl0')
   type(ctrl_out),save :: o_SWupTOA      = ctrl_out((/ 1, 4, 10, 10, 10 /),'SWupTOA')
@@ -108,4 +115,6 @@
   type(ctrl_out),save :: o_SWdnTOA      = ctrl_out((/ 1, 4, 10, 10, 10 /),'SWdnTOA')
   type(ctrl_out),save :: o_SWdnTOAclr   = ctrl_out((/ 1, 4, 10, 10, 10 /),'SWdnTOAclr')
+  type(ctrl_out),save :: o_nettop       = ctrl_out((/ 1, 4, 10, 10, 10 /),'nettop')
+
   type(ctrl_out),save :: o_SWup200      = ctrl_out((/ 1, 10, 10, 10, 10 /),'SWup200')
   type(ctrl_out),save :: o_SWup200clr   = ctrl_out((/ 10, 1, 10, 10, 10 /),'SWup200clr')
@@ -123,49 +132,51 @@
   type(ctrl_out),save :: o_LWdn200      = ctrl_out((/ 1, 10, 10, 10, 10 /),'LWdn200')
   type(ctrl_out),save :: o_LWdn200clr   = ctrl_out((/ 1, 10, 10, 10, 10 /),'LWdn200clr')
-  type(ctrl_out),save :: o_sols         = ctrl_out((/ 1, 1, 10, 1, 10 /),'sols')
+  type(ctrl_out),save :: o_sols         = ctrl_out((/ 1, 1, 10, 10, 10 /),'sols')
   type(ctrl_out),save :: o_sols0        = ctrl_out((/ 1, 5, 10, 10, 10 /),'sols0')
-  type(ctrl_out),save :: o_soll         = ctrl_out((/ 1, 1, 10, 1, 10 /),'soll')
+  type(ctrl_out),save :: o_soll         = ctrl_out((/ 1, 1, 10, 10, 10 /),'soll')
   type(ctrl_out),save :: o_soll0        = ctrl_out((/ 1, 5, 10, 10, 10 /),'soll0')
   type(ctrl_out),save :: o_radsol       = ctrl_out((/ 1, 1, 10, 10, 10 /),'radsol')
-  type(ctrl_out),save :: o_SWupSFC      = ctrl_out((/ 1, 4, 10, 10, 10 /),'SWupSFC')
+  type(ctrl_out),save :: o_SWupSFC      = ctrl_out((/ 1, 4, 10, 5, 10 /),'SWupSFC')
   type(ctrl_out),save :: o_SWupSFCclr   = ctrl_out((/ 1, 4, 10, 10, 10 /),'SWupSFCclr')
-  type(ctrl_out),save :: o_SWdnSFC      = ctrl_out((/ 1, 1, 10, 10, 10 /),'SWdnSFC') 
-  type(ctrl_out),save :: o_SWdnSFCclr   = ctrl_out((/ 1, 4, 10, 10, 10 /),'SWdnSFCclr')
+  type(ctrl_out),save :: o_SWdnSFC      = ctrl_out((/ 1, 1, 10, 5, 10 /),'SWdnSFC') 
+  type(ctrl_out),save :: o_SWdnSFCclr   = ctrl_out((/ 1, 4, 10, 5, 10 /),'SWdnSFCclr')
   type(ctrl_out),save :: o_LWupSFC      = ctrl_out((/ 1, 4, 10, 10, 10 /),'LWupSFC')
-  type(ctrl_out),save :: o_LWupSFCclr   = ctrl_out((/ 1, 4, 10, 10, 10 /),'LWupSFCclr')
-  type(ctrl_out),save :: o_LWdnSFC      = ctrl_out((/ 1, 4, 10, 10, 10 /),'LWdnSFC')
-  type(ctrl_out),save :: o_LWdnSFCclr   = ctrl_out((/ 1, 4, 10, 10, 10 /),'LWdnSFCclr')
-  type(ctrl_out),save :: o_bils         = ctrl_out((/ 1, 2, 10, 1, 10 /),'bils')
-  type(ctrl_out),save :: o_sens         = ctrl_out((/ 1, 1, 10, 1, 1 /),'sens')
-  type(ctrl_out),save :: o_fder         = ctrl_out((/ 1, 2, 10, 1, 10 /),'fder')
+  type(ctrl_out),save :: o_LWupSFCclr   = ctrl_out((/ 1, 4, 10, 5, 10 /),'LWupSFCclr')
+  type(ctrl_out),save :: o_LWdnSFC      = ctrl_out((/ 1, 4, 10, 5, 10 /),'LWdnSFC')
+  type(ctrl_out),save :: o_LWdnSFCclr   = ctrl_out((/ 1, 4, 10, 5, 10 /),'LWdnSFCclr')
+  type(ctrl_out),save :: o_bils         = ctrl_out((/ 1, 2, 10, 5, 10 /),'bils')
+  type(ctrl_out),save :: o_sens         = ctrl_out((/ 1, 1, 10, 5, 10 /),'sens')
+  type(ctrl_out),save :: o_fder         = ctrl_out((/ 1, 2, 10, 10, 10 /),'fder')
   type(ctrl_out),save :: o_ffonte       = ctrl_out((/ 1, 10, 10, 10, 10 /),'ffonte')
   type(ctrl_out),save :: o_fqcalving    = ctrl_out((/ 1, 10, 10, 10, 10 /),'fqcalving')
   type(ctrl_out),save :: o_fqfonte      = ctrl_out((/ 1, 10, 10, 10, 10 /),'fqfonte')
 
-  type(ctrl_out),save,dimension(4) :: o_taux_srf     = (/ ctrl_out((/ 1, 4, 10, 1, 10 /),'taux_ter'), &
-                                                 ctrl_out((/ 1, 4, 10, 1, 10 /),'taux_lic'), &
-                                                 ctrl_out((/ 1, 4, 10, 1, 10 /),'taux_oce'), &
-                                                 ctrl_out((/ 1, 4, 10, 1, 10 /),'taux_sic') /)
-
-  type(ctrl_out),save,dimension(4) :: o_tauy_srf     = (/ ctrl_out((/ 1, 4, 10, 1, 10 /),'tauy_ter'), &
-                                                 ctrl_out((/ 1, 4, 10, 1, 10 /),'tauy_lic'), &
-                                                 ctrl_out((/ 1, 4, 10, 1, 10 /),'tauy_oce'), &
-                                                 ctrl_out((/ 1, 4, 10, 1, 10 /),'tauy_sic') /)
-
-
-  type(ctrl_out),save,dimension(4) :: o_pourc_srf    = (/ ctrl_out((/ 1, 4, 10, 1, 10 /),'pourc_ter'), &
-                                                 ctrl_out((/ 1, 4, 10, 1, 10 /),'pourc_lic'), &
-                                                 ctrl_out((/ 1, 4, 10, 1, 10 /),'pourc_oce'), &
-                                                 ctrl_out((/ 1, 4, 10, 1, 10 /),'pourc_sic') /)     
-
-  type(ctrl_out),save,dimension(4) :: o_fract_srf    = (/ ctrl_out((/ 1, 4, 10, 1, 10 /),'fract_ter'), &
-                                                 ctrl_out((/ 1, 4, 10, 1, 10 /),'fract_lic'), &
-                                                 ctrl_out((/ 1, 4, 10, 1, 10 /),'fract_oce'), &
-                                                 ctrl_out((/ 1, 4, 10, 1, 10 /),'fract_sic') /)
-
-  type(ctrl_out),save,dimension(4) :: o_tsol_srf     = (/ ctrl_out((/ 1, 4, 10, 1, 10 /),'tsol_ter'), &
-                                                 ctrl_out((/ 1, 4, 10, 1, 10 /),'tsol_lic'), &
-                                                 ctrl_out((/ 1, 4, 10, 1, 10 /),'tsol_oce'), &
-                                                 ctrl_out((/ 1, 4, 10, 1, 10 /),'tsol_sic') /)
+  type(ctrl_out),save :: o_taux         = ctrl_out((/ 1, 10, 10, 10, 10 /),'taux')
+  type(ctrl_out),save :: o_tauy         = ctrl_out((/ 1, 10, 10, 10, 10 /),'tauy')
+  type(ctrl_out),save,dimension(4) :: o_taux_srf     = (/ ctrl_out((/ 1, 4, 10, 10, 10 /),'taux_ter'), &
+                                                 ctrl_out((/ 1, 4, 10, 10, 10 /),'taux_lic'), &
+                                                 ctrl_out((/ 1, 4, 10, 10, 10 /),'taux_oce'), &
+                                                 ctrl_out((/ 1, 4, 10, 10, 10 /),'taux_sic') /)
+
+  type(ctrl_out),save,dimension(4) :: o_tauy_srf     = (/ ctrl_out((/ 1, 4, 10, 10, 10 /),'tauy_ter'), &
+                                                 ctrl_out((/ 1, 4, 10, 10, 10 /),'tauy_lic'), &
+                                                 ctrl_out((/ 1, 4, 10, 10, 10 /),'tauy_oce'), &
+                                                 ctrl_out((/ 1, 4, 10, 10, 10 /),'tauy_sic') /)
+
+
+  type(ctrl_out),save,dimension(4) :: o_pourc_srf    = (/ ctrl_out((/ 1, 4, 10, 10, 10 /),'pourc_ter'), &
+                                                 ctrl_out((/ 1, 4, 10, 10, 10 /),'pourc_lic'), &
+                                                 ctrl_out((/ 1, 4, 10, 10, 10 /),'pourc_oce'), &
+                                                 ctrl_out((/ 1, 4, 10, 10, 10 /),'pourc_sic') /)     
+
+  type(ctrl_out),save,dimension(4) :: o_fract_srf    = (/ ctrl_out((/ 1, 4, 10, 10, 10 /),'fract_ter'), &
+                                                 ctrl_out((/ 1, 4, 10, 10, 10 /),'fract_lic'), &
+                                                 ctrl_out((/ 1, 4, 10, 10, 10 /),'fract_oce'), &
+                                                 ctrl_out((/ 1, 4, 10, 10, 10 /),'fract_sic') /)
+
+  type(ctrl_out),save,dimension(4) :: o_tsol_srf     = (/ ctrl_out((/ 1, 4, 10, 10, 10 /),'tsol_ter'), &
+                                                 ctrl_out((/ 1, 4, 10, 10, 10 /),'tsol_lic'), &
+                                                 ctrl_out((/ 1, 4, 10, 10, 10 /),'tsol_oce'), &
+                                                 ctrl_out((/ 1, 4, 10, 10, 10 /),'tsol_sic') /)
 
   type(ctrl_out),save,dimension(4) :: o_sens_srf     = (/ ctrl_out((/ 1, 4, 10, 1, 10 /),'sens_ter'), &
@@ -200,10 +211,10 @@
 
 
-  type(ctrl_out),save :: o_cdrm         = ctrl_out((/ 1, 10, 10, 1, 10 /),'cdrm')
+  type(ctrl_out),save :: o_cdrm         = ctrl_out((/ 1, 10, 10, 10, 10 /),'cdrm')
   type(ctrl_out),save :: o_cdrh         = ctrl_out((/ 1, 10, 10, 1, 10 /),'cdrh')
   type(ctrl_out),save :: o_cldl         = ctrl_out((/ 1, 1, 10, 10, 10 /),'cldl')
   type(ctrl_out),save :: o_cldm         = ctrl_out((/ 1, 1, 10, 10, 10 /),'cldm')
   type(ctrl_out),save :: o_cldh         = ctrl_out((/ 1, 1, 10, 10, 10 /),'cldh')
-  type(ctrl_out),save :: o_cldt         = ctrl_out((/ 1, 1, 2, 10, 10 /),'cldt')
+  type(ctrl_out),save :: o_cldt         = ctrl_out((/ 1, 1, 2, 5, 10 /),'cldt')
   type(ctrl_out),save :: o_cldq         = ctrl_out((/ 1, 1, 10, 10, 10 /),'cldq')
   type(ctrl_out),save :: o_lwp          = ctrl_out((/ 1, 5, 10, 10, 10 /),'lwp')
@@ -215,31 +226,31 @@
  
   type(ctrl_out),save :: o_cape         = ctrl_out((/ 1, 10, 10, 10, 10 /),'cape')
-  type(ctrl_out),save :: o_pbase        = ctrl_out((/ 1, 10, 10, 10, 10 /),'pbase')
-  type(ctrl_out),save :: o_ptop         = ctrl_out((/ 1, 4, 10, 10, 10 /),'ptop')
+  type(ctrl_out),save :: o_pbase        = ctrl_out((/ 1, 5, 10, 10, 10 /),'pbase')
+  type(ctrl_out),save :: o_ptop         = ctrl_out((/ 1, 5, 10, 10, 10 /),'ptop')
   type(ctrl_out),save :: o_fbase        = ctrl_out((/ 1, 10, 10, 10, 10 /),'fbase')
   type(ctrl_out),save :: o_prw          = ctrl_out((/ 1, 1, 10, 10, 10 /),'prw')
 
-  type(ctrl_out),save :: o_s_pblh       = ctrl_out((/ 1, 10, 10, 1, 1 /),'s_pblh')
-  type(ctrl_out),save :: o_s_pblt       = ctrl_out((/ 1, 10, 10, 1, 1 /),'s_pblt')
-  type(ctrl_out),save :: o_s_lcl        = ctrl_out((/ 1, 10, 10, 1, 10 /),'s_lcl')
-  type(ctrl_out),save :: o_s_capCL      = ctrl_out((/ 1, 10, 10, 1, 10 /),'s_capCL')
-  type(ctrl_out),save :: o_s_oliqCL     = ctrl_out((/ 1, 10, 10, 1, 10 /),'s_oliqCL')
-  type(ctrl_out),save :: o_s_cteiCL     = ctrl_out((/ 1, 10, 10, 1, 1 /),'s_cteiCL')
-  type(ctrl_out),save :: o_s_therm      = ctrl_out((/ 1, 10, 10, 1, 1 /),'s_therm')
-  type(ctrl_out),save :: o_s_trmb1      = ctrl_out((/ 1, 10, 10, 1, 10 /),'s_trmb1')
-  type(ctrl_out),save :: o_s_trmb2      = ctrl_out((/ 1, 10, 10, 1, 10 /),'s_trmb2')
-  type(ctrl_out),save :: o_s_trmb3      = ctrl_out((/ 1, 10, 10, 1, 10 /),'s_trmb3')
+  type(ctrl_out),save :: o_s_pblh       = ctrl_out((/ 1, 10, 10, 10, 1 /),'s_pblh')
+  type(ctrl_out),save :: o_s_pblt       = ctrl_out((/ 1, 10, 10, 10, 1 /),'s_pblt')
+  type(ctrl_out),save :: o_s_lcl        = ctrl_out((/ 1, 10, 10, 10, 10 /),'s_lcl')
+  type(ctrl_out),save :: o_s_capCL      = ctrl_out((/ 1, 10, 10, 10, 10 /),'s_capCL')
+  type(ctrl_out),save :: o_s_oliqCL     = ctrl_out((/ 1, 10, 10, 10, 10 /),'s_oliqCL')
+  type(ctrl_out),save :: o_s_cteiCL     = ctrl_out((/ 1, 10, 10, 10, 1 /),'s_cteiCL')
+  type(ctrl_out),save :: o_s_therm      = ctrl_out((/ 1, 10, 10, 10, 1 /),'s_therm')
+  type(ctrl_out),save :: o_s_trmb1      = ctrl_out((/ 1, 10, 10, 10, 10 /),'s_trmb1')
+  type(ctrl_out),save :: o_s_trmb2      = ctrl_out((/ 1, 10, 10, 10, 10 /),'s_trmb2')
+  type(ctrl_out),save :: o_s_trmb3      = ctrl_out((/ 1, 10, 10, 10, 10 /),'s_trmb3')
 
   type(ctrl_out),save :: o_slab_bils    = ctrl_out((/ 1, 1, 10, 10, 10 /),'slab_bils_oce')
 
-  type(ctrl_out),save :: o_ale_bl       = ctrl_out((/ 1, 1, 1, 1, 10 /),'ale_bl')
-  type(ctrl_out),save :: o_alp_bl       = ctrl_out((/ 1, 1, 1, 1, 10 /),'alp_bl')
-  type(ctrl_out),save :: o_ale_wk       = ctrl_out((/ 1, 1, 1, 1, 10 /),'ale_wk')
-  type(ctrl_out),save :: o_alp_wk       = ctrl_out((/ 1, 1, 1, 1, 10 /),'alp_wk')
-
-  type(ctrl_out),save :: o_ale          = ctrl_out((/ 1, 1, 1, 1, 10 /),'ale')
-  type(ctrl_out),save :: o_alp          = ctrl_out((/ 1, 1, 1, 1, 10 /),'alp')
-  type(ctrl_out),save :: o_cin          = ctrl_out((/ 1, 1, 1, 1, 10 /),'cin')
-  type(ctrl_out),save :: o_wape         = ctrl_out((/ 1, 1, 1, 1, 10 /),'wape')
+  type(ctrl_out),save :: o_ale_bl       = ctrl_out((/ 1, 1, 1, 10, 10 /),'ale_bl')
+  type(ctrl_out),save :: o_alp_bl       = ctrl_out((/ 1, 1, 1, 10, 10 /),'alp_bl')
+  type(ctrl_out),save :: o_ale_wk       = ctrl_out((/ 1, 1, 1, 10, 10 /),'ale_wk')
+  type(ctrl_out),save :: o_alp_wk       = ctrl_out((/ 1, 1, 1, 10, 10 /),'alp_wk')
+
+  type(ctrl_out),save :: o_ale          = ctrl_out((/ 1, 1, 1, 10, 10 /),'ale')
+  type(ctrl_out),save :: o_alp          = ctrl_out((/ 1, 1, 1, 10, 10 /),'alp')
+  type(ctrl_out),save :: o_cin          = ctrl_out((/ 1, 1, 1, 10, 10 /),'cin')
+  type(ctrl_out),save :: o_wape         = ctrl_out((/ 1, 1, 1, 10, 10 /),'wape')
 
 
@@ -306,5 +317,7 @@
   type(ctrl_out),save :: o_cldtau       = ctrl_out((/ 10, 5, 10, 10, 10 /),'cldtau')                     
   type(ctrl_out),save :: o_cldemi       = ctrl_out((/ 10, 5, 10, 10, 10 /),'cldemi')
-  type(ctrl_out),save :: o_rh2m         = ctrl_out((/ 10, 5, 10, 10, 10 /),'rh2m')
+  type(ctrl_out),save :: o_rh2m         = ctrl_out((/ 5, 5, 10, 10, 10 /),'rh2m')
+  type(ctrl_out),save :: o_rh2m_min     = ctrl_out((/ 10, 5, 10, 10, 10 /),'rh2m_min')
+  type(ctrl_out),save :: o_rh2m_max     = ctrl_out((/ 10, 5, 10, 10, 10 /),'rh2m_max')
   type(ctrl_out),save :: o_qsat2m       = ctrl_out((/ 10, 5, 10, 10, 10 /),'qsat2m')
   type(ctrl_out),save :: o_tpot         = ctrl_out((/ 10, 5, 10, 10, 10 /),'tpot')
@@ -331,11 +344,11 @@
   type(ctrl_out),save :: o_snowl        = ctrl_out((/ 10, 1, 10, 10, 10 /),'snowl')
   type(ctrl_out),save :: o_cape_max     = ctrl_out((/ 10, 1, 10, 10, 10 /),'cape_max')
-  type(ctrl_out),save :: o_solldown     = ctrl_out((/ 10, 1, 10, 1, 10 /),'solldown')
-
-  type(ctrl_out),save :: o_dtsvdfo      = ctrl_out((/ 10, 10, 10, 1, 10 /),'dtsvdfo')
-  type(ctrl_out),save :: o_dtsvdft      = ctrl_out((/ 10, 10, 10, 1, 10 /),'dtsvdft')
-  type(ctrl_out),save :: o_dtsvdfg      = ctrl_out((/ 10, 10, 10, 1, 10 /),'dtsvdfg')
-  type(ctrl_out),save :: o_dtsvdfi      = ctrl_out((/ 10, 10, 10, 1, 10 /),'dtsvdfi')
-  type(ctrl_out),save :: o_rugs         = ctrl_out((/ 10, 10, 10, 1, 1 /),'rugs')
+  type(ctrl_out),save :: o_solldown     = ctrl_out((/ 10, 1, 10, 10, 10 /),'solldown')
+
+  type(ctrl_out),save :: o_dtsvdfo      = ctrl_out((/ 10, 10, 10, 10, 10 /),'dtsvdfo')
+  type(ctrl_out),save :: o_dtsvdft      = ctrl_out((/ 10, 10, 10, 10, 10 /),'dtsvdft')
+  type(ctrl_out),save :: o_dtsvdfg      = ctrl_out((/ 10, 10, 10, 10, 10 /),'dtsvdfg')
+  type(ctrl_out),save :: o_dtsvdfi      = ctrl_out((/ 10, 10, 10, 10, 10 /),'dtsvdfi')
+  type(ctrl_out),save :: o_rugs         = ctrl_out((/ 10, 10, 10, 10, 10 /),'rugs')
 
   type(ctrl_out),save :: o_topswad      = ctrl_out((/ 2, 10, 10, 10, 10 /),'topswad')
@@ -344,49 +357,77 @@
   type(ctrl_out),save :: o_solswai      = ctrl_out((/ 2, 10, 10, 10, 10 /),'solswai')
 
-  type(ctrl_out),save,dimension(10) :: o_tausumaero  = (/ ctrl_out((/ 4, 10, 10, 10, 10 /),'OD550_ASBCM'), &
-                                                     ctrl_out((/ 4, 10, 10, 10, 10 /),'OD550_ASPOMM'), &
-                                                     ctrl_out((/ 4, 10, 10, 10, 10 /),'OD550_ASSO4M'), &
-                                                     ctrl_out((/ 4, 10, 10, 10, 10 /),'OD550_CSSO4M'), &
-                                                     ctrl_out((/ 4, 10, 10, 10, 10 /),'OD550_SSSSM'), &
-                                                     ctrl_out((/ 4, 10, 10, 10, 10 /),'OD550_ASSSM'), &
-                                                     ctrl_out((/ 4, 10, 10, 10, 10 /),'OD550_CSSSM'), &
-                                                     ctrl_out((/ 4, 10, 10, 10, 10 /),'OD550_CIDUSTM'), &
-                                                     ctrl_out((/ 4, 10, 10, 10, 10 /),'OD550_AIBCM'), &
-                                                     ctrl_out((/ 4, 10, 10, 10, 10 /),'OD550_AIPOMM') /)
-
-  type(ctrl_out),save :: o_swtoaas_nat      = ctrl_out((/ 4, 10, 10, 10, 10 /),'swtoaas_nat')
-  type(ctrl_out),save :: o_swsrfas_nat      = ctrl_out((/ 4, 10, 10, 10, 10 /),'swsrfas_nat')
-  type(ctrl_out),save :: o_swtoacs_nat      = ctrl_out((/ 4, 10, 10, 10, 10 /),'swtoacs_nat')
-  type(ctrl_out),save :: o_swsrfcs_nat      = ctrl_out((/ 4, 10, 10, 10, 10 /),'swsrfcs_nat')
-
-  type(ctrl_out),save :: o_swtoaas_ant      = ctrl_out((/ 4, 10, 10, 10, 10 /),'swtoaas_ant')
-  type(ctrl_out),save :: o_swsrfas_ant      = ctrl_out((/ 4, 10, 10, 10, 10 /),'swsrfas_ant')
-  type(ctrl_out),save :: o_swtoacs_ant      = ctrl_out((/ 4, 10, 10, 10, 10 /),'swtoacs_ant')
-  type(ctrl_out),save :: o_swsrfcs_ant      = ctrl_out((/ 4, 10, 10, 10, 10 /),'swsrfcs_ant')
-
-  type(ctrl_out),save :: o_swtoacf_nat      = ctrl_out((/ 4, 10, 10, 10, 10 /),'swtoacf_nat')
-  type(ctrl_out),save :: o_swsrfcf_nat      = ctrl_out((/ 4, 10, 10, 10, 10 /),'swsrfcf_nat')
-  type(ctrl_out),save :: o_swtoacf_ant      = ctrl_out((/ 4, 10, 10, 10, 10 /),'swtoacf_ant')
-  type(ctrl_out),save :: o_swsrfcf_ant      = ctrl_out((/ 4, 10, 10, 10, 10 /),'swsrfcf_ant')
-  type(ctrl_out),save :: o_swtoacf_zero      = ctrl_out((/ 4, 10, 10, 10, 10 /),'swtoacf_zero')
-  type(ctrl_out),save :: o_swsrfcf_zero      = ctrl_out((/ 4, 10, 10, 10, 10 /),'swsrfcf_zero')
+  type(ctrl_out),save,dimension(10) :: o_tausumaero  = (/ ctrl_out((/ 4, 4, 10, 10, 10 /),'OD550_ASBCM'), &
+                                                     ctrl_out((/ 4, 4, 10, 10, 10 /),'OD550_ASPOMM'), &
+                                                     ctrl_out((/ 4, 4, 10, 10, 10 /),'OD550_ASSO4M'), &
+                                                     ctrl_out((/ 4, 4, 10, 10, 10 /),'OD550_CSSO4M'), &
+                                                     ctrl_out((/ 4, 4, 10, 10, 10 /),'OD550_SSSSM'), &
+                                                     ctrl_out((/ 4, 4, 10, 10, 10 /),'OD550_ASSSM'), &
+                                                     ctrl_out((/ 4, 4, 10, 10, 10 /),'OD550_CSSSM'), &
+                                                     ctrl_out((/ 4, 4, 10, 10, 10 /),'OD550_CIDUSTM'), &
+                                                     ctrl_out((/ 4, 4, 10, 10, 10 /),'OD550_AIBCM'), &
+                                                     ctrl_out((/ 4, 4, 10, 10, 10 /),'OD550_AIPOMM') /)
+
+  type(ctrl_out),save :: o_od550aer         = ctrl_out((/ 4, 4, 10, 10, 10 /),'od550aer')
+  type(ctrl_out),save :: o_od865aer         = ctrl_out((/ 4, 4, 10, 10, 10 /),'od865aer')
+  type(ctrl_out),save :: o_absvisaer        = ctrl_out((/ 4, 4, 10, 10, 10 /),'absvisaer')
+  type(ctrl_out),save :: o_od550lt1aer      = ctrl_out((/ 4, 4, 10, 10, 10 /),'od550lt1aer')
+
+  type(ctrl_out),save :: o_sconcso4         = ctrl_out((/ 4, 4, 10, 10, 10 /),'sconcso4')
+  type(ctrl_out),save :: o_sconcoa          = ctrl_out((/ 4, 4, 10, 10, 10 /),'sconcoa')
+  type(ctrl_out),save :: o_sconcbc          = ctrl_out((/ 4, 4, 10, 10, 10 /),'sconcbc')
+  type(ctrl_out),save :: o_sconcss          = ctrl_out((/ 4, 4, 10, 10, 10 /),'sconcss')
+  type(ctrl_out),save :: o_sconcdust        = ctrl_out((/ 4, 4, 10, 10, 10 /),'sconcdust')
+  type(ctrl_out),save :: o_concso4          = ctrl_out((/ 4, 4, 10, 10, 10 /),'concso4')
+  type(ctrl_out),save :: o_concoa           = ctrl_out((/ 4, 4, 10, 10, 10 /),'concoa')
+  type(ctrl_out),save :: o_concbc           = ctrl_out((/ 4, 4, 10, 10, 10 /),'concbc')
+  type(ctrl_out),save :: o_concss           = ctrl_out((/ 4, 4, 10, 10, 10 /),'concss')
+  type(ctrl_out),save :: o_concdust         = ctrl_out((/ 4, 4, 10, 10, 10 /),'concdust')
+  type(ctrl_out),save :: o_loadso4          = ctrl_out((/ 4, 4, 10, 10, 10 /),'loadso4')
+  type(ctrl_out),save :: o_loadoa           = ctrl_out((/ 4, 4, 10, 10, 10 /),'loadoa')
+  type(ctrl_out),save :: o_loadbc           = ctrl_out((/ 4, 4, 10, 10, 10 /),'loadbc')
+  type(ctrl_out),save :: o_loadss           = ctrl_out((/ 4, 4, 10, 10, 10 /),'loadss')
+  type(ctrl_out),save :: o_loaddust         = ctrl_out((/ 4, 4, 10, 10, 10 /),'loaddust')
+
+  type(ctrl_out),save :: o_swtoaas_nat      = ctrl_out((/ 4, 4, 10, 10, 10 /),'swtoaas_nat')
+  type(ctrl_out),save :: o_swsrfas_nat      = ctrl_out((/ 4, 4, 10, 10, 10 /),'swsrfas_nat')
+  type(ctrl_out),save :: o_swtoacs_nat      = ctrl_out((/ 4, 4, 10, 10, 10 /),'swtoacs_nat')
+  type(ctrl_out),save :: o_swsrfcs_nat      = ctrl_out((/ 4, 4, 10, 10, 10 /),'swsrfcs_nat')
+
+  type(ctrl_out),save :: o_swtoaas_ant      = ctrl_out((/ 4, 4, 10, 10, 10 /),'swtoaas_ant')
+  type(ctrl_out),save :: o_swsrfas_ant      = ctrl_out((/ 4, 4, 10, 10, 10 /),'swsrfas_ant')
+  type(ctrl_out),save :: o_swtoacs_ant      = ctrl_out((/ 4, 4, 10, 10, 10 /),'swtoacs_ant')
+  type(ctrl_out),save :: o_swsrfcs_ant      = ctrl_out((/ 4, 4, 10, 10, 10 /),'swsrfcs_ant')
+
+  type(ctrl_out),save :: o_swtoacf_nat      = ctrl_out((/ 4, 4, 10, 10, 10 /),'swtoacf_nat')
+  type(ctrl_out),save :: o_swsrfcf_nat      = ctrl_out((/ 4, 4, 10, 10, 10 /),'swsrfcf_nat')
+  type(ctrl_out),save :: o_swtoacf_ant      = ctrl_out((/ 4, 4, 10, 10, 10 /),'swtoacf_ant')
+  type(ctrl_out),save :: o_swsrfcf_ant      = ctrl_out((/ 4, 4, 10, 10, 10 /),'swsrfcf_ant')
+  type(ctrl_out),save :: o_swtoacf_zero     = ctrl_out((/ 4, 4, 10, 10, 10 /),'swtoacf_zero')
+  type(ctrl_out),save :: o_swsrfcf_zero     = ctrl_out((/ 4, 4, 10, 10, 10 /),'swsrfcf_zero')
+
+  type(ctrl_out),save :: o_cldncl          = ctrl_out((/ 4, 4, 10, 10, 10 /),'cldncl')
+  type(ctrl_out),save :: o_reffclwtop      = ctrl_out((/ 4, 4, 10, 10, 10 /),'reffclwtop')
+  type(ctrl_out),save :: o_cldnvi          = ctrl_out((/ 4, 4, 10, 10, 10 /),'cldnvi')
+  type(ctrl_out),save :: o_lcc             = ctrl_out((/ 4, 4, 10, 10, 10 /),'lcc')
 
 
 !!!!!!!!!!!!!!!!!!!!!! 3D !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
+  type(ctrl_out),save :: o_ec550aer     = ctrl_out((/ 4,  4, 10, 10, 1 /),'ec550aer')
   type(ctrl_out),save :: o_lwcon        = ctrl_out((/ 2, 5, 10, 10, 1 /),'lwcon')
   type(ctrl_out),save :: o_iwcon        = ctrl_out((/ 2, 5, 10, 10, 10 /),'iwcon')
-  type(ctrl_out),save :: o_temp         = ctrl_out((/ 2, 3, 4, 1, 1 /),'temp')
-  type(ctrl_out),save :: o_theta        = ctrl_out((/ 2, 3, 4, 1, 1 /),'theta')
-  type(ctrl_out),save :: o_ovap         = ctrl_out((/ 2, 3, 4, 1, 1 /),'ovap')
-  type(ctrl_out),save :: o_ovapinit         = ctrl_out((/ 2, 3, 4, 1, 1 /),'ovapinit')
+  type(ctrl_out),save :: o_temp         = ctrl_out((/ 2, 3, 4, 10, 1 /),'temp')
+  type(ctrl_out),save :: o_theta        = ctrl_out((/ 2, 3, 4, 10, 1 /),'theta')
+  type(ctrl_out),save :: o_ovap         = ctrl_out((/ 2, 3, 4, 10, 1 /),'ovap')
+  type(ctrl_out),save :: o_ovapinit         = ctrl_out((/ 2, 3, 10, 10, 1 /),'ovapinit')
   type(ctrl_out),save :: o_wvapp        = ctrl_out((/ 2, 10, 10, 10, 10 /),'wvapp')
-  type(ctrl_out),save :: o_geop         = ctrl_out((/ 2, 3, 10, 1, 1 /),'geop')
-  type(ctrl_out),save :: o_vitu         = ctrl_out((/ 2, 3, 4, 1, 1 /),'vitu')
-  type(ctrl_out),save :: o_vitv         = ctrl_out((/ 2, 3, 4, 1, 1 /),'vitv')
-  type(ctrl_out),save :: o_vitw         = ctrl_out((/ 2, 3, 10, 10, 1 /),'vitw')
-  type(ctrl_out),save :: o_pres         = ctrl_out((/ 2, 3, 10, 1, 1 /),'pres')
+  type(ctrl_out),save :: o_geop         = ctrl_out((/ 2, 3, 10, 10, 1 /),'geop')
+  type(ctrl_out),save :: o_vitu         = ctrl_out((/ 2, 3, 4, 5, 1 /),'vitu')
+  type(ctrl_out),save :: o_vitv         = ctrl_out((/ 2, 3, 4, 5, 1 /),'vitv')
+  type(ctrl_out),save :: o_vitw         = ctrl_out((/ 2, 3, 10, 5, 1 /),'vitw')
+  type(ctrl_out),save :: o_pres         = ctrl_out((/ 2, 3, 10, 10, 1 /),'pres')
+  type(ctrl_out),save :: o_paprs        = ctrl_out((/ 2, 3, 10, 10, 1 /),'paprs')
   type(ctrl_out),save :: o_rneb         = ctrl_out((/ 2, 5, 10, 10, 1 /),'rneb')
   type(ctrl_out),save :: o_rnebcon      = ctrl_out((/ 2, 5, 10, 10, 1 /),'rnebcon')
-  type(ctrl_out),save :: o_rhum         = ctrl_out((/ 2, 10, 10, 10, 10 /),'rhum')
+  type(ctrl_out),save :: o_rhum         = ctrl_out((/ 2, 5, 10, 10, 10 /),'rhum')
   type(ctrl_out),save :: o_ozone        = ctrl_out((/ 2, 10, 10, 10, 10 /),'ozone')
   type(ctrl_out),save :: o_ozone_light        = ctrl_out((/ 2, 10, 10, 10, 10 /),'ozone_daylight')
@@ -400,4 +441,10 @@
   type(ctrl_out),save :: o_re           =ctrl_out((/ 5, 10, 10, 10, 10 /),'re')
   type(ctrl_out),save :: o_fl           =ctrl_out((/ 5, 10, 10, 10, 10 /),'fl')
+  type(ctrl_out),save :: o_scdnc        =ctrl_out((/ 4,  4, 10, 10, 1 /),'scdnc')
+  type(ctrl_out),save :: o_reffclws     =ctrl_out((/ 4,  4, 10, 10, 1 /),'reffclws')
+  type(ctrl_out),save :: o_reffclwc     =ctrl_out((/ 4,  4, 10, 10, 1 /),'reffclwc')
+  type(ctrl_out),save :: o_lcc3d        =ctrl_out((/ 4,  4, 10, 10, 1 /),'lcc3d')
+  type(ctrl_out),save :: o_lcc3dcon     =ctrl_out((/ 4,  4, 10, 10, 1 /),'lcc3dcon')
+  type(ctrl_out),save :: o_lcc3dstra    =ctrl_out((/ 4,  4, 10, 10, 1 /),'lcc3dstra')
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
@@ -412,11 +459,11 @@
                                                      ctrl_out((/ 3, 10, 10, 10, 10 /),'ages_sic') /)
 
-  type(ctrl_out),save,dimension(4) :: o_rugs_srf     = (/ ctrl_out((/ 3, 4, 10, 1, 10 /),'rugs_ter'), &
-                                                     ctrl_out((/ 3, 4, 10, 1, 10 /),'rugs_lic'), &
-                                                     ctrl_out((/ 3, 4, 10, 1, 10 /),'rugs_oce'), &
-                                                     ctrl_out((/ 3, 4, 10, 1, 10 /),'rugs_sic') /)
-
-  type(ctrl_out),save :: o_albs         = ctrl_out((/ 3, 10, 10, 1, 10 /),'albs')
-  type(ctrl_out),save :: o_albslw       = ctrl_out((/ 3, 10, 10, 1, 10 /),'albslw')
+  type(ctrl_out),save,dimension(4) :: o_rugs_srf     = (/ ctrl_out((/ 3, 4, 10, 10, 10 /),'rugs_ter'), &
+                                                     ctrl_out((/ 3, 4, 10, 10, 10 /),'rugs_lic'), &
+                                                     ctrl_out((/ 3, 4, 10, 10, 10 /),'rugs_oce'), &
+                                                     ctrl_out((/ 3, 4, 10, 10, 10 /),'rugs_sic') /)
+
+  type(ctrl_out),save :: o_albs         = ctrl_out((/ 3, 10, 10, 10, 10 /),'albs')
+  type(ctrl_out),save :: o_albslw       = ctrl_out((/ 3, 10, 10, 10, 10 /),'albslw')
 
   type(ctrl_out),save :: o_clwcon       = ctrl_out((/ 4, 10, 10, 10, 10 /),'clwcon')
@@ -424,4 +471,6 @@
   type(ctrl_out),save :: o_dnwd         = ctrl_out((/ 4, 10, 10, 10, 10 /),'dnwd')
   type(ctrl_out),save :: o_dnwd0        = ctrl_out((/ 4, 10, 10, 10, 10 /),'dnwd0')
+  type(ctrl_out),save :: o_mc           = ctrl_out((/ 4, 10, 10, 10, 10 /),'mc')
+  type(ctrl_out),save :: o_ftime_con    = ctrl_out((/ 4, 10, 10, 10, 10 /),'ftime_con')
   type(ctrl_out),save :: o_dtdyn        = ctrl_out((/ 4, 10, 10, 10, 1 /),'dtdyn')
   type(ctrl_out),save :: o_dqdyn        = ctrl_out((/ 4, 10, 10, 10, 1 /),'dqdyn')
@@ -444,6 +493,6 @@
   type(ctrl_out),save :: o_dtlschr      = ctrl_out((/ 4, 10, 10, 10, 10 /),'dtlschr')
   type(ctrl_out),save :: o_dqlsc        = ctrl_out((/ 4, 10, 10, 10, 10 /),'dqlsc')
-  type(ctrl_out),save :: o_dtvdf        = ctrl_out((/ 4, 10, 10, 1, 10 /),'dtvdf')
-  type(ctrl_out),save :: o_dqvdf        = ctrl_out((/ 4, 10, 10, 1, 10 /),'dqvdf')
+  type(ctrl_out),save :: o_dtvdf        = ctrl_out((/ 4, 10, 10, 10, 10 /),'dtvdf')
+  type(ctrl_out),save :: o_dqvdf        = ctrl_out((/ 4, 10, 10, 10, 10 /),'dqvdf')
   type(ctrl_out),save :: o_dteva        = ctrl_out((/ 4, 10, 10, 10, 10 /),'dteva')
   type(ctrl_out),save :: o_dqeva        = ctrl_out((/ 4, 10, 10, 10, 10 /),'dqeva')
@@ -455,4 +504,5 @@
   type(ctrl_out),save :: o_w_th         = ctrl_out((/ 4, 10, 10, 10, 10 /),'w_th')
   type(ctrl_out),save :: o_lambda_th    = ctrl_out((/ 10, 10, 10, 10, 10 /),'lambda_th')
+  type(ctrl_out),save :: o_ftime_th     = ctrl_out((/ 10, 10, 10, 10, 10 /),'ftime_th')
   type(ctrl_out),save :: o_q_th         = ctrl_out((/ 4, 10, 10, 10, 10 /),'q_th')
   type(ctrl_out),save :: o_a_th         = ctrl_out((/ 4, 10, 10, 10, 10 /),'a_th')
@@ -772,4 +822,6 @@
  CALL histdef2d(iff,o_pluc%flag,o_pluc%name, "Convective Precip.", "kg/(s*m2)")
  CALL histdef2d(iff,o_snow%flag,o_snow%name, "Snow fall", "kg/(s*m2)" )
+ CALL histdef2d(iff,o_msnow%flag,o_msnow%name, "Surface snow amount", "kg/m2" )
+ CALL histdef2d(iff,o_fsnow%flag,o_fsnow%name, "Surface snow area fraction", "-" )
  CALL histdef2d(iff,o_evap%flag,o_evap%name, "Evaporat", "kg/(s*m2)" ) 
  CALL histdef2d(iff,o_tops%flag,o_tops%name, "Solar rad. at TOA", "W/m2") 
@@ -781,4 +833,5 @@
  CALL histdef2d(iff,o_SWdnTOA%flag,o_SWdnTOA%name, "SWdn at TOA", "W/m2" )
  CALL histdef2d(iff,o_SWdnTOAclr%flag,o_SWdnTOAclr%name, "SWdn clear sky at TOA", "W/m2") 
+ CALL histdef2d(iff,o_nettop%flag,o_nettop%name, "Net dn radiatif flux at TOA", "W/m2")
  CALL histdef2d(iff,o_SWup200%flag,o_SWup200%name, "SWup at 200mb", "W/m2" ) 
  CALL histdef2d(iff,o_SWup200clr%flag,o_SWup200clr%name, "SWup clear sky at 200mb", "W/m2")
@@ -809,4 +862,7 @@
  CALL histdef2d(iff,o_fqfonte%flag,o_fqfonte%name, "Land ice melt", "kg/m2/s") 
 
+ CALL histdef2d(iff,o_taux%flag,o_taux%name, "Zonal wind stress","Pa")
+ CALL histdef2d(iff,o_tauy%flag,o_tauy%name, "Meridional wind stress","Pa")
+
      DO nsrf = 1, nbsrf
  CALL histdef2d(iff,o_pourc_srf(nsrf)%flag,o_pourc_srf(nsrf)%name,"% "//clnsurf(nsrf),"%")
@@ -816,4 +872,5 @@
  CALL histdef2d(iff,o_tsol_srf(nsrf)%flag,o_tsol_srf(nsrf)%name,"Temperature "//clnsurf(nsrf),"K")
  CALL histdef2d(iff,o_u10m_srf(nsrf)%flag,o_u10m_srf(nsrf)%name,"Vent Zonal 10m "//clnsurf(nsrf),"m/s")
+ CALL histdef2d(iff,o_evap_srf(nsrf)%flag,o_evap_srf(nsrf)%name,"evaporation at surface "//clnsurf(nsrf),"kg/(s*m2)")
  CALL histdef2d(iff,o_v10m_srf(nsrf)%flag,o_v10m_srf(nsrf)%name,"Vent meredien 10m "//clnsurf(nsrf),"m/s")
  CALL histdef2d(iff,o_t2m_srf(nsrf)%flag,o_t2m_srf(nsrf)%name,"Temp 2m "//clnsurf(nsrf),"K")
@@ -840,4 +897,27 @@
 
 IF (new_aod .AND. (.NOT. aerosol_couple)) THEN
+
+  CALL histdef2d(iff,o_od550aer%flag,o_od550aer%name, "Total aerosol optical depth at 550nm", "-")
+  CALL histdef2d(iff,o_od865aer%flag,o_od865aer%name, "Total aerosol optical depth at 870nm", "-")
+  CALL histdef2d(iff,o_absvisaer%flag,o_absvisaer%name, "Absorption aerosol visible optical depth", "-")
+  CALL histdef2d(iff,o_od550lt1aer%flag,o_od550lt1aer%name, "Fine mode optical depth", "-")
+
+
+  CALL histdef2d(iff,o_sconcso4%flag,o_sconcso4%name,"Surface Concentration of Sulfate ","kg/m3")
+  CALL histdef2d(iff,o_sconcoa%flag,o_sconcoa%name,"Surface Concentration of Organic Aerosol ","kg/m3")
+  CALL histdef2d(iff,o_sconcbc%flag,o_sconcbc%name,"Surface Concentration of Black Carbon ","kg/m3")
+  CALL histdef2d(iff,o_sconcss%flag,o_sconcss%name,"Surface Concentration of Sea Salt ","kg/m3")
+  CALL histdef2d(iff,o_sconcdust%flag,o_sconcdust%name,"Surface Concentration of Dust ","kg/m3")
+  CALL histdef3d(iff,o_concso4%flag,o_concso4%name,"Concentration of Sulfate ","kg/m3")
+  CALL histdef3d(iff,o_concoa%flag,o_concoa%name,"Concentration of Organic Aerosol ","kg/m3")
+  CALL histdef3d(iff,o_concbc%flag,o_concbc%name,"Concentration of Black Carbon ","kg/m3")
+  CALL histdef3d(iff,o_concss%flag,o_concss%name,"Concentration of Sea Salt ","kg/m3")
+  CALL histdef3d(iff,o_concdust%flag,o_concdust%name,"Concentration of Dust ","kg/m3")
+  CALL histdef2d(iff,o_loadso4%flag,o_loadso4%name,"Column Load of Sulfate ","kg/m2")
+  CALL histdef2d(iff,o_loadoa%flag,o_loadoa%name,"Column Load of Organic Aerosol ","kg/m2")
+  CALL histdef2d(iff,o_loadbc%flag,o_loadbc%name,"Column Load of Black Carbon ","kg/m2")
+  CALL histdef2d(iff,o_loadss%flag,o_loadss%name,"Column Load of Sea Salt ","kg/m2")
+  CALL histdef2d(iff,o_loaddust%flag,o_loaddust%name,"Column Load of Dust ","kg/m2")
+
   DO naero = 1, naero_spc
   CALL histdef2d(iff,o_tausumaero(naero)%flag,o_tausumaero(naero)%name,"Aerosol Optical depth at 550 nm "//name_aero(naero),"1")
@@ -873,4 +953,15 @@
   CALL histdef2d(iff,o_topswai%flag,o_topswai%name, "AIE at TOA", "W/m2")
   CALL histdef2d(iff,o_solswai%flag,o_solswai%name, "AIE at SFR", "W/m2")
+!Cloud droplet number concentration
+  CALL histdef3d(iff,o_scdnc%flag,o_scdnc%name, "Cloud droplet number concentration","m-3")
+  CALL histdef2d(iff,o_cldncl%flag,o_cldncl%name, "CDNC at top of liquid water cloud", "m-3")
+  CALL histdef3d(iff,o_reffclws%flag,o_reffclws%name, "Stratiform Cloud Droplet Effective Radius","m")
+  CALL histdef3d(iff,o_reffclwc%flag,o_reffclwc%name, "Convective Cloud Droplet Effective Radius","m")
+  CALL histdef2d(iff,o_cldnvi%flag,o_cldnvi%name, "Column Integrated Cloud Droplet Number", "m-2")
+  CALL histdef3d(iff,o_lcc3d%flag,o_lcc3d%name, "Cloud liquid fraction","1")
+  CALL histdef3d(iff,o_lcc3dcon%flag,o_lcc3dcon%name, "Convective cloud liquid fraction","1")
+  CALL histdef3d(iff,o_lcc3dstra%flag,o_lcc3dstra%name, "Stratiform cloud liquid fraction","1")
+  CALL histdef2d(iff,o_lcc%flag,o_lcc%name, "Cloud liquid fraction at top of cloud","1")
+  CALL histdef2d(iff,o_reffclwtop%flag,o_reffclwtop%name, "Droplet effective radius at top of liquid water cloud", "m")
  ENDIF
 
@@ -909,4 +1000,12 @@
  CALL histdef3d(iff,o_dnwd%flag,o_dnwd%name, "saturated downdraft", "kg/m2/s")
  CALL histdef3d(iff,o_dnwd0%flag,o_dnwd0%name, "unsat. downdraft", "kg/m2/s")
+ CALL histdef3d(iff,o_mc%flag,o_mc%name, "Convective mass flux", "kg/m2/s")
+   type_ecri(1) = 'inst(X)'
+   type_ecri(2) = 'inst(X)'
+   type_ecri(3) = 'inst(X)'
+   type_ecri(4) = 'inst(X)'
+   type_ecri(5) = 'inst(X)'
+ CALL histdef2d(iff,o_ftime_con%flag,o_ftime_con%name, "Fraction of time convection Occurs", " ")
+   type_ecri(:) = type_ecri_files(:)
      ENDIF !iflag_con .GE. 3
 
@@ -965,4 +1064,17 @@
  CALL histdef2d(iff,o_dthmin%flag,o_dthmin%name, "dTheta mini", "K/m")
  CALL histdef2d(iff,o_rh2m%flag,o_rh2m%name, "Relative humidity at 2m", "%" )
+   type_ecri(1) = 't_min(X)'
+   type_ecri(2) = 't_min(X)'
+   type_ecri(3) = 't_min(X)'
+   type_ecri(4) = 't_min(X)'
+   type_ecri(5) = 't_min(X)'
+ CALL histdef2d(iff,o_rh2m_min%flag,o_rh2m_min%name, "Min Relative humidity at 2m", "%" )
+   type_ecri(1) = 't_max(X)'
+   type_ecri(2) = 't_max(X)'
+   type_ecri(3) = 't_max(X)'
+   type_ecri(4) = 't_max(X)'
+   type_ecri(5) = 't_max(X)'
+ CALL histdef2d(iff,o_rh2m_max%flag,o_rh2m_max%name, "Max Relative humidity at 2m", "%" )
+   type_ecri(:) = type_ecri_files(:)
  CALL histdef2d(iff,o_qsat2m%flag,o_qsat2m%name, "Saturant humidity at 2m", "%")
  CALL histdef2d(iff,o_tpot%flag,o_tpot%name, "Surface air potential temperature", "K")
@@ -981,4 +1093,5 @@
 
 ! Champs 3D:
+ CALL histdef3d(iff,o_ec550aer%flag,o_ec550aer%name, "Extinction at 550nm", "m^-1")
  CALL histdef3d(iff,o_lwcon%flag,o_lwcon%name, "Cloud liquid water content", "kg/kg")
  CALL histdef3d(iff,o_iwcon%flag,o_iwcon%name, "Cloud ice water content", "kg/kg")
@@ -992,4 +1105,5 @@
  CALL histdef3d(iff,o_vitw%flag,o_vitw%name, "Vertical wind", "Pa/s" )
  CALL histdef3d(iff,o_pres%flag,o_pres%name, "Air pressure", "Pa" )
+ CALL histdef3d(iff,o_paprs%flag,o_paprs%name, "Air pressure Inter-Couches", "Pa" )
  CALL histdef3d(iff,o_rneb%flag,o_rneb%name, "Cloud fraction", "-")
  CALL histdef3d(iff,o_rnebcon%flag,o_rnebcon%name, "Convective Cloud Fraction", "-")
@@ -1078,4 +1192,5 @@
  CALL histdef3d(iff,o_w_th%flag,o_w_th%name,"Thermal plume vertical velocity","m/s")
  CALL histdef3d(iff,o_lambda_th%flag,o_lambda_th%name,"Thermal plume vertical velocity","m/s")
+ CALL histdef2d(iff,o_ftime_th%flag,o_ftime_th%name,"Fraction of time Shallow convection occurs"," ")
  CALL histdef3d(iff,o_q_th%flag,o_q_th%name, "Thermal plume total humidity", "kg/kg")
  CALL histdef3d(iff,o_a_th%flag,o_a_th%name, "Thermal plume fraction", "")
Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phys_output_write.h
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phys_output_write.h	(revision 1346)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phys_output_write.h	(revision 1347)
@@ -146,4 +146,12 @@
        ENDIF
 
+       IF (o_msnow%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),o_msnow%name,itau_w,snow_o)
+       ENDIF
+
+       IF (o_fsnow%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),o_fsnow%name,itau_w,zfra_o)
+       ENDIF
+
        IF (o_evap%flag(iff)<=lev_files(iff)) THEN
       CALL histwrite_phy(nid_files(iff),o_evap%name,itau_w,evap)
@@ -190,4 +198,10 @@
        ENDIF
 
+       IF (o_nettop%flag(iff)<=lev_files(iff)) THEN
+      zx_tmp_fi2d(:) = topsw(:)-toplw(:)
+      CALL histwrite_phy(nid_files(iff),
+     $                  o_nettop%name,itau_w,zx_tmp_fi2d)
+       ENDIF
+
        IF (o_SWup200%flag(iff)<=lev_files(iff)) THEN
       CALL histwrite_phy(nid_files(iff),o_SWup200%name,itau_w,SWup200)
@@ -319,4 +333,22 @@
        CALL histwrite_phy(nid_files(iff),
      $                   o_fqfonte%name,itau_w,zxfqfonte)
+       ENDIF
+
+       IF (o_taux%flag(iff)<=lev_files(iff)) THEN
+         zx_tmp_fi2d=0.
+         do nsrf=1,nbsrf
+          zx_tmp_fi2d(:)=zx_tmp_fi2d(:)+pctsrf(:,nsrf)*fluxu(:,1,nsrf)
+         enddo
+         CALL histwrite_phy(nid_files(iff),
+     $                   o_taux%name,itau_w,zx_tmp_fi2d)
+       ENDIF
+
+       IF (o_tauy%flag(iff)<=lev_files(iff)) THEN
+         zx_tmp_fi2d=0.
+         do nsrf=1,nbsrf
+          zx_tmp_fi2d(:)=zx_tmp_fi2d(:)+pctsrf(:,nsrf)*fluxv(:,1,nsrf)
+         enddo
+         CALL histwrite_phy(nid_files(iff),
+     $                   o_tauy%name,itau_w,zx_tmp_fi2d)
        ENDIF
 
@@ -377,4 +409,10 @@
       ENDIF
 
+      IF (o_evap_srf(nsrf)%flag(iff)<=lev_files(iff)) THEN
+      zx_tmp_fi2d(1 : klon) = fevap(1 : klon, nsrf)
+      CALL histwrite_phy(nid_files(iff),o_evap_srf(nsrf)%name,
+     $           itau_w,zx_tmp_fi2d)
+      ENDIF
+
        IF (o_sens_srf(nsrf)%flag(iff)<=lev_files(iff)) THEN
        zx_tmp_fi2d(1 : klon) = fluxt( 1 : klon, 1, nsrf)
@@ -532,4 +570,19 @@
        ENDIF
 
+       IF (o_ftime_con%flag(iff)<=lev_files(iff)) THEN
+        zx_tmp_fi2d=float(itau_con)/float(itap)
+      CALL histwrite_phy(nid_files(iff),o_ftime_con%name,
+     s                   itau_w,zx_tmp_fi2d)
+       ENDIF
+
+       IF (o_mc%flag(iff)<=lev_files(iff)) THEN
+        if(iflag_thermals.gt.1)then
+         zx_tmp_fi3d=dnwd+dnwd0+upwd+fm_therm
+        else
+         zx_tmp_fi3d=dnwd+dnwd0+upwd
+        endif 
+      CALL histwrite_phy(nid_files(iff),o_mc%name,itau_w,zx_tmp_fi3d)
+       ENDIF
+      
       ENDIF !iflag_con .GE. 3
 
@@ -780,4 +833,21 @@
       ENDIF
 
+      IF (o_rh2m_min%flag(iff)<=lev_files(iff)) THEN
+      DO i=1, klon
+       zx_tmp_fi2d(i)=MIN(100.,rh2m(i)*100.)
+      ENDDO
+      CALL histwrite_phy(nid_files(iff),o_rh2m_min%name,
+     s               itau_w,zx_tmp_fi2d)
+      ENDIF
+
+      IF (o_rh2m_max%flag(iff)<=lev_files(iff)) THEN
+      DO i=1, klon
+       zx_tmp_fi2d(i)=MIN(100.,rh2m(i)*100.)
+      ENDDO
+      CALL histwrite_phy(nid_files(iff),o_rh2m_max%name,
+     s              itau_w,zx_tmp_fi2d)
+      ENDIF
+
+
       IF (o_qsat2m%flag(iff)<=lev_files(iff)) THEN
       CALL histwrite_phy(nid_files(iff),o_qsat2m%name,itau_w,qsat2m)
@@ -1127,4 +1197,9 @@
         IF (o_pres%flag(iff)<=lev_files(iff)) THEN
       CALL histwrite_phy(nid_files(iff),o_pres%name,itau_w,pplay)
+        ENDIF
+
+        IF (o_paprs%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),o_paprs%name,
+     s                    itau_w,paprs(:,1:klev))
         ENDIF
 
@@ -1307,4 +1382,11 @@
 
        IF (iflag_thermals.gt.1) THEN
+        IF (o_ftime_th%flag(iff)<=lev_files(iff)) THEN
+! Pour l instant 0 a y reflichir pour les thermiques
+         zx_tmp_fi2d=0. 
+        CALL histwrite_phy(nid_files(iff),o_ftime_th%name,
+     s                     itau_w,zx_tmp_fi2d)
+        ENDIF
+
         IF (o_f_th%flag(iff)<=lev_files(iff)) THEN
         CALL histwrite_phy(nid_files(iff),o_f_th%name,itau_w,fm_therm)
@@ -1336,5 +1418,5 @@
        ENDIF
 
-       ENDIF !iflag_thermals
+      ENDIF !iflag_thermals
 
        IF (o_f0_th%flag(iff)<=lev_files(iff)) THEN
Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F	(revision 1346)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F	(revision 1347)
@@ -29,4 +29,5 @@
       USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
       USE phys_state_var_mod ! Variables sauvegardees de la physique
+      USE phys_output_var_mod ! Variables pour les ecritures des sorties
       USE fonte_neige_mod, ONLY  : fonte_neige_get_vars
       USE phys_output_mod
Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/readaerosol_interp.F90
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/readaerosol_interp.F90	(revision 1346)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/readaerosol_interp.F90	(revision 1347)
@@ -1,5 +1,5 @@
 ! $Id$
 !
-SUBROUTINE readaerosol_interp(id_aero, itap, pdtphys, r_day, first, pplay, paprs, t_seri, mass_out, pi_mass_out)
+SUBROUTINE readaerosol_interp(id_aero, itap, pdtphys, r_day, first, pplay, paprs, t_seri, mass_out, pi_mass_out, load_src)
 !
 ! This routine will return the mass concentration at actual day(mass_out) and 
@@ -46,4 +46,5 @@
   REAL, INTENT(OUT) :: mass_out(klon,klev)    ! Mass of aerosol (monthly mean data,from file) [ug AIBCM/m3]
   REAL, INTENT(OUT) :: pi_mass_out(klon,klev) ! Mass of preindustrial aerosol (monthly mean data,from file) [ug AIBCM/m3]
+  REAL, INTENT(OUT) :: load_src(klon) ! Load of aerosol (monthly mean data,from file) [kg/m3]
 !      
 ! Local Variables:
@@ -61,5 +62,5 @@
   REAL                              :: volm      ! Volyme de melange [kg/kg]
   REAL, DIMENSION(klon)             :: psurf_day, pi_psurf_day
-  REAL, DIMENSION(klon)             :: load_src, pi_load_src  ! Mass load at source grid
+  REAL, DIMENSION(klon)             :: pi_load_src  ! Mass load at source grid
   REAL, DIMENSION(klon)             :: load_tgt, load_tgt_test
   REAL, DIMENSION(klon,klev)        :: delp ! pressure difference in each model layer
Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/readaerosol_optic.F90
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/readaerosol_optic.F90	(revision 1346)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/readaerosol_optic.F90	(revision 1347)
@@ -14,4 +14,7 @@
   USE dimphy
   USE aero_mod
+  USE phys_local_var_mod, only: sconcso4,sconcoa,sconcbc,sconcss,sconcdust, &
+      concso4,concoa,concbc,concss,concdust,loadso4,loadoa,loadbc,loadss,loaddust, &
+      load_tmp1,load_tmp2,load_tmp3,load_tmp4,load_tmp5,load_tmp6,load_tmp7
   IMPLICIT NONE
 
@@ -76,7 +79,8 @@
        flag_aerosol .EQ. 6 ) THEN 
 
-     CALL readaerosol_interp(id_ASSO4M, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sulfate, sulfate_pi)
+     CALL readaerosol_interp(id_ASSO4M, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sulfate, sulfate_pi,loadso4)
   ELSE
      sulfate(:,:) = 0. ; sulfate_pi(:,:) = 0.
+     loadso4=0.
   END IF
 
@@ -86,9 +90,11 @@
 
      ! Get bc aerosol distribution 
-     CALL readaerosol_interp(id_ASBCM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, bcsol, bcsol_pi )
-     CALL readaerosol_interp(id_AIBCM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, bcins, bcins_pi )
+     CALL readaerosol_interp(id_ASBCM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, bcsol, bcsol_pi, load_tmp1 )
+     CALL readaerosol_interp(id_AIBCM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, bcins, bcins_pi, load_tmp2 )
+     loadbc(:)=load_tmp1(:)+load_tmp2(:)
   ELSE
      bcsol(:,:) = 0. ; bcsol_pi(:,:) = 0.
      bcins(:,:) = 0. ; bcins_pi(:,:) = 0.
+     loadbc=0.
   END IF
 
@@ -98,9 +104,11 @@
        flag_aerosol .EQ. 6 ) THEN
 
-     CALL readaerosol_interp(id_ASPOMM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomsol, pomsol_pi)
-     CALL readaerosol_interp(id_AIPOMM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomins, pomins_pi)
+     CALL readaerosol_interp(id_ASPOMM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomsol, pomsol_pi, load_tmp3)
+     CALL readaerosol_interp(id_AIPOMM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, pomins, pomins_pi, load_tmp4)
+     loadoa(:)=load_tmp3(:)+load_tmp4(:)
   ELSE
      pomsol(:,:) = 0. ; pomsol_pi(:,:) = 0.
      pomins(:,:) = 0. ; pomins_pi(:,:) = 0.
+     loadoa=0.
   END IF
 
@@ -110,12 +118,13 @@
       flag_aerosol .EQ. 6 ) THEN 
 
-      CALL readaerosol_interp(id_SSSSM ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sssupco, sssupco_pi) 
-      CALL readaerosol_interp(id_CSSSM ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sscoarse,sscoarse_pi) 
-      CALL readaerosol_interp(id_ASSSM ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, ssacu, ssacu_pi) 
-
+      CALL readaerosol_interp(id_SSSSM ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sssupco, sssupco_pi, load_tmp5) 
+      CALL readaerosol_interp(id_CSSSM ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, sscoarse,sscoarse_pi, load_tmp6) 
+      CALL readaerosol_interp(id_ASSSM ,itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, ssacu, ssacu_pi, load_tmp7) 
+     loadss(:)=load_tmp5(:)+load_tmp6(:)+load_tmp7(:)
   ELSE
      sscoarse(:,:) = 0. ; sscoarse_pi(:,:) = 0. 
      ssacu(:,:)    = 0. ; ssacu_pi(:,:) = 0. 
      sssupco(:,:)  = 0. ; sssupco_pi = 0. 
+     loadss=0.
   ENDIF
 
@@ -124,8 +133,9 @@
       flag_aerosol .EQ. 6 ) THEN 
 
-      CALL readaerosol_interp(id_CIDUSTM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, cidust, cidust_pi) 
+      CALL readaerosol_interp(id_CIDUSTM, itap, pdtphys, rjourvrai, debut, pplay, paprs, t_seri, cidust, cidust_pi, loaddust) 
 
   ELSE
       cidust(:,:) = 0. ; cidust_pi(:,:) = 0. 
+      loaddust=0.
   ENDIF
 
@@ -198,3 +208,17 @@
   END IF
 
+
+! Diagnostics calculation for CMIP5 protocol
+  sconcso4(:)=m_allaer(:,1,id_ASSO4M)*1.e-9
+  sconcoa(:)=(m_allaer(:,1,id_ASPOMM)+m_allaer(:,1,id_AIPOMM))*1.e-9
+  sconcbc(:)=(m_allaer(:,1,id_ASBCM)+m_allaer(:,1,id_AIBCM))*1.e-9
+  sconcss(:)=(m_allaer(:,1,id_ASSSM)+m_allaer(:,1,id_CSSSM)+m_allaer(:,1,id_SSSSM))*1.e-9
+  sconcdust(:)=m_allaer(:,1,id_CIDUSTM)*1.e-9
+  concso4(:,:)=m_allaer(:,:,id_ASSO4M)*1.e-9
+  concoa(:,:)=(m_allaer(:,:,id_ASPOMM)+m_allaer(:,:,id_AIPOMM))*1.e-9
+  concbc(:,:)=(m_allaer(:,:,id_ASBCM)+m_allaer(:,:,id_AIBCM))*1.e-9
+  concss(:,:)=(m_allaer(:,:,id_ASSSM)+m_allaer(:,:,id_CSSSM)+m_allaer(:,:,id_SSSSM))*1.e-9
+  concdust(:,:)=m_allaer(:,:,id_CIDUSTM)*1.e-9
+
+
 END SUBROUTINE readaerosol_optic
