Index: LMDZ5/trunk/libf/phylmd/cosp/mrgrnk.F90.new
===================================================================
--- LMDZ5/trunk/libf/phylmd/cosp/mrgrnk.F90.new	(revision 2432)
+++ 	(revision )
@@ -1,611 +1,0 @@
-! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $
-! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/quickbeam/mrgrnk.f90 $
-Module m_mrgrnk
-Integer, Parameter :: kdp = selected_real_kind(15)
-public :: mrgrnk
-private :: kdp
-private :: R_mrgrnk, I_mrgrnk, D_mrgrnk
-interface mrgrnk
-  module procedure D_mrgrnk, R_mrgrnk, I_mrgrnk
-end interface mrgrnk
-contains
-
-Subroutine D_mrgrnk (XDONT, IRNGT)
-! __________________________________________________________
-!   MRGRNK = Merge-sort ranking of an array
-!   For performance reasons, the first 2 passes are taken
-!   out of the standard loop, and use dedicated coding.
-! __________________________________________________________
-! __________________________________________________________
-      Real (kind=kdp), Dimension (:), Intent (In) :: XDONT
-      Integer*4, Dimension (:), Intent (Out) :: IRNGT
-! __________________________________________________________
-      Real (kind=kdp) :: XVALA, XVALB
-!
-      Integer, Dimension (SIZE(IRNGT)) :: JWRKT
-      Integer :: LMTNA, LMTNC, IRNG1, IRNG2
-      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
-!
-      NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
-      Select Case (NVAL)
-      Case (:0)
-         Return
-      Case (1)
-         IRNGT (1) = 1
-         Return
-      Case Default
-         Continue
-      End Select
-!
-!  Fill-in the index array, creating ordered couples
-!
-      Do IIND = 2, NVAL, 2
-         If (XDONT(IIND-1) <= XDONT(IIND)) Then
-            IRNGT (IIND-1) = IIND - 1
-            IRNGT (IIND) = IIND
-         Else
-            IRNGT (IIND-1) = IIND
-            IRNGT (IIND) = IIND - 1
-         End If
-      End Do
-      If (Modulo(NVAL, 2) /= 0) Then
-         IRNGT (NVAL) = NVAL
-      End If
-!
-!  We will now have ordered subsets A - B - A - B - ...
-!  and merge A and B couples into     C   -   C   - ...
-!
-      LMTNA = 2
-      LMTNC = 4
-!
-!  First iteration. The length of the ordered subsets goes from 2 to 4
-!
-      Do
-         If (NVAL <= 2) Exit
-!
-!   Loop on merges of A and B into C
-!
-         Do IWRKD = 0, NVAL - 1, 4
-            If ((IWRKD+4) > NVAL) Then
-               If ((IWRKD+2) >= NVAL) Exit
-!
-!   1 2 3
-!
-               If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
-!
-!   1 3 2
-!
-               If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
-                  IRNG2 = IRNGT (IWRKD+2)
-                  IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
-                  IRNGT (IWRKD+3) = IRNG2
-!
-!   3 1 2
-!
-               Else
-                  IRNG1 = IRNGT (IWRKD+1)
-                  IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
-                  IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
-                  IRNGT (IWRKD+2) = IRNG1
-               End If
-               Exit
-            End If
-!
-!   1 2 3 4
-!
-            If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
-!
-!   1 3 x x
-!
-            If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
-               IRNG2 = IRNGT (IWRKD+2)
-               IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
-               If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
-!   1 3 2 4
-                  IRNGT (IWRKD+3) = IRNG2
-               Else
-!   1 3 4 2
-                  IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
-                  IRNGT (IWRKD+4) = IRNG2
-               End If
-!
-!   3 x x x
-!
-            Else
-               IRNG1 = IRNGT (IWRKD+1)
-               IRNG2 = IRNGT (IWRKD+2)
-               IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
-               If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
-                  IRNGT (IWRKD+2) = IRNG1
-                  If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
-!   3 1 2 4
-                     IRNGT (IWRKD+3) = IRNG2
-                  Else
-!   3 1 4 2
-                     IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
-                     IRNGT (IWRKD+4) = IRNG2
-                  End If
-               Else
-!   3 4 1 2
-                  IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
-                  IRNGT (IWRKD+3) = IRNG1
-                  IRNGT (IWRKD+4) = IRNG2
-               End If
-            End If
-         End Do
-!
-!  The Cs become As and Bs
-!
-         LMTNA = 4
-         Exit
-      End Do
-!
-!  Iteration loop. Each time, the length of the ordered subsets
-!  is doubled.
-!
-      Do
-         If (LMTNA >= NVAL) Exit
-         IWRKF = 0
-         LMTNC = 2 * LMTNC
-!
-!   Loop on merges of A and B into C
-!
-         Do
-            IWRK = IWRKF
-            IWRKD = IWRKF + 1
-            JINDA = IWRKF + LMTNA
-            IWRKF = IWRKF + LMTNC
-            If (IWRKF >= NVAL) Then
-               If (JINDA >= NVAL) Exit
-               IWRKF = NVAL
-            End If
-            IINDA = 1
-            IINDB = JINDA + 1
-!
-!   Shortcut for the case when the max of A is smaller
-!   than the min of B. This line may be activated when the
-!   initial set is already close to sorted.
-!
-!          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
-!
-!  One steps in the C subset, that we build in the final rank array
-!
-!  Make a copy of the rank array for the merge iteration
-!
-            JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
-!
-            XVALA = XDONT (JWRKT(IINDA))
-            XVALB = XDONT (IRNGT(IINDB))
-!
-            Do
-               IWRK = IWRK + 1
-!
-!  We still have unprocessed values in both A and B
-!
-               If (XVALA > XVALB) Then
-                  IRNGT (IWRK) = IRNGT (IINDB)
-                  IINDB = IINDB + 1
-                  If (IINDB > IWRKF) Then
-!  Only A still with unprocessed values
-                     IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
-                     Exit
-                  End If
-                  XVALB = XDONT (IRNGT(IINDB))
-               Else
-                  IRNGT (IWRK) = JWRKT (IINDA)
-                  IINDA = IINDA + 1
-                  If (IINDA > LMTNA) Exit! Only B still with unprocessed values
-                  XVALA = XDONT (JWRKT(IINDA))
-               End If
-!
-            End Do
-         End Do
-!
-!  The Cs become As and Bs
-!
-         LMTNA = 2 * LMTNA
-      End Do
-!
-      Return
-!
-End Subroutine D_mrgrnk
-
-Subroutine R_mrgrnk (XDONT, IRNGT)
-! __________________________________________________________
-!   MRGRNK = Merge-sort ranking of an array
-!   For performance reasons, the first 2 passes are taken
-!   out of the standard loop, and use dedicated coding.
-! __________________________________________________________
-! _________________________________________________________
-      Real, Dimension (:), Intent (In) :: XDONT
-      Integer*4, Dimension (:), Intent (Out) :: IRNGT
-! __________________________________________________________
-      Real :: XVALA, XVALB
-!
-      Integer, Dimension (SIZE(IRNGT)) :: JWRKT
-      Integer :: LMTNA, LMTNC, IRNG1, IRNG2
-      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
-!
-      NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
-      Select Case (NVAL)
-      Case (:0)
-         Return
-      Case (1)
-         IRNGT (1) = 1
-         Return
-      Case Default
-         Continue
-      End Select
-!
-!  Fill-in the index array, creating ordered couples
-!
-      Do IIND = 2, NVAL, 2
-         If (XDONT(IIND-1) <= XDONT(IIND)) Then
-            IRNGT (IIND-1) = IIND - 1
-            IRNGT (IIND) = IIND
-         Else
-            IRNGT (IIND-1) = IIND
-            IRNGT (IIND) = IIND - 1
-         End If
-      End Do
-      If (Modulo(NVAL, 2) /= 0) Then
-         IRNGT (NVAL) = NVAL
-      End If
-!
-!  We will now have ordered subsets A - B - A - B - ...
-!  and merge A and B couples into     C   -   C   - ...
-!
-      LMTNA = 2
-      LMTNC = 4
-!
-!  First iteration. The length of the ordered subsets goes from 2 to 4
-!
-      Do
-         If (NVAL <= 2) Exit
-!
-!   Loop on merges of A and B into C
-!
-         Do IWRKD = 0, NVAL - 1, 4
-            If ((IWRKD+4) > NVAL) Then
-               If ((IWRKD+2) >= NVAL) Exit
-!
-!   1 2 3
-!
-               If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
-!
-!   1 3 2
-!
-               If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
-                  IRNG2 = IRNGT (IWRKD+2)
-                  IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
-                  IRNGT (IWRKD+3) = IRNG2
-!
-!   3 1 2
-!
-               Else
-                  IRNG1 = IRNGT (IWRKD+1)
-                  IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
-                  IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
-                  IRNGT (IWRKD+2) = IRNG1
-               End If
-               Exit
-            End If
-!
-!   1 2 3 4
-!
-            If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
-!
-!   1 3 x x
-!
-            If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
-               IRNG2 = IRNGT (IWRKD+2)
-               IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
-               If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
-!   1 3 2 4
-                  IRNGT (IWRKD+3) = IRNG2
-               Else
-!   1 3 4 2
-                  IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
-                  IRNGT (IWRKD+4) = IRNG2
-               End If
-!
-!   3 x x x
-!
-            Else
-               IRNG1 = IRNGT (IWRKD+1)
-               IRNG2 = IRNGT (IWRKD+2)
-               IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
-               If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
-                  IRNGT (IWRKD+2) = IRNG1
-                  If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
-!   3 1 2 4
-                     IRNGT (IWRKD+3) = IRNG2
-                  Else
-!   3 1 4 2
-                     IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
-                     IRNGT (IWRKD+4) = IRNG2
-                  End If
-               Else
-!   3 4 1 2
-                  IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
-                  IRNGT (IWRKD+3) = IRNG1
-                  IRNGT (IWRKD+4) = IRNG2
-               End If
-            End If
-         End Do
-!
-!  The Cs become As and Bs
-!
-         LMTNA = 4
-         Exit
-      End Do
-!
-!  Iteration loop. Each time, the length of the ordered subsets
-!  is doubled.
-!
-      Do
-         If (LMTNA >= NVAL) Exit
-         IWRKF = 0
-         LMTNC = 2 * LMTNC
-!
-!   Loop on merges of A and B into C
-!
-         Do
-            IWRK = IWRKF
-            IWRKD = IWRKF + 1
-            JINDA = IWRKF + LMTNA
-            IWRKF = IWRKF + LMTNC
-            If (IWRKF >= NVAL) Then
-               If (JINDA >= NVAL) Exit
-               IWRKF = NVAL
-            End If
-            IINDA = 1
-            IINDB = JINDA + 1
-!
-!   Shortcut for the case when the max of A is smaller
-!   than the min of B. This line may be activated when the
-!   initial set is already close to sorted.
-!
-!          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
-!
-!  One steps in the C subset, that we build in the final rank array
-!
-!  Make a copy of the rank array for the merge iteration
-!
-            JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
-!
-            XVALA = XDONT (JWRKT(IINDA))
-            XVALB = XDONT (IRNGT(IINDB))
-!
-            Do
-               IWRK = IWRK + 1
-!
-!  We still have unprocessed values in both A and B
-!
-               If (XVALA > XVALB) Then
-                  IRNGT (IWRK) = IRNGT (IINDB)
-                  IINDB = IINDB + 1
-                  If (IINDB > IWRKF) Then
-!  Only A still with unprocessed values
-                     IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
-                     Exit
-                  End If
-                  XVALB = XDONT (IRNGT(IINDB))
-               Else
-                  IRNGT (IWRK) = JWRKT (IINDA)
-                  IINDA = IINDA + 1
-                  If (IINDA > LMTNA) Exit! Only B still with unprocessed values
-                  XVALA = XDONT (JWRKT(IINDA))
-               End If
-!
-            End Do
-         End Do
-!
-!  The Cs become As and Bs
-!
-         LMTNA = 2 * LMTNA
-      End Do
-!
-      Return
-!
-End Subroutine R_mrgrnk
-Subroutine I_mrgrnk (XDONT, IRNGT)
-! __________________________________________________________
-!   MRGRNK = Merge-sort ranking of an array
-!   For performance reasons, the first 2 passes are taken
-!   out of the standard loop, and use dedicated coding.
-! __________________________________________________________
-! __________________________________________________________
-      Integer, Dimension (:), Intent (In)  :: XDONT
-      Integer*4, Dimension (:), Intent (Out) :: IRNGT
-! __________________________________________________________
-      Integer :: XVALA, XVALB
-!
-      Integer, Dimension (SIZE(IRNGT)) :: JWRKT
-      Integer :: LMTNA, LMTNC, IRNG1, IRNG2
-      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
-!
-      NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
-      Select Case (NVAL)
-      Case (:0)
-         Return
-      Case (1)
-         IRNGT (1) = 1
-         Return
-      Case Default
-         Continue
-      End Select
-!
-!  Fill-in the index array, creating ordered couples
-!
-      Do IIND = 2, NVAL, 2
-         If (XDONT(IIND-1) <= XDONT(IIND)) Then
-            IRNGT (IIND-1) = IIND - 1
-            IRNGT (IIND) = IIND
-         Else
-            IRNGT (IIND-1) = IIND
-            IRNGT (IIND) = IIND - 1
-         End If
-      End Do
-      If (Modulo(NVAL, 2) /= 0) Then
-         IRNGT (NVAL) = NVAL
-      End If
-!
-!  We will now have ordered subsets A - B - A - B - ...
-!  and merge A and B couples into     C   -   C   - ...
-!
-      LMTNA = 2
-      LMTNC = 4
-!
-!  First iteration. The length of the ordered subsets goes from 2 to 4
-!
-      Do
-         If (NVAL <= 2) Exit
-!
-!   Loop on merges of A and B into C
-!
-         Do IWRKD = 0, NVAL - 1, 4
-            If ((IWRKD+4) > NVAL) Then
-               If ((IWRKD+2) >= NVAL) Exit
-!
-!   1 2 3
-!
-               If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
-!
-!   1 3 2
-!
-               If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
-                  IRNG2 = IRNGT (IWRKD+2)
-                  IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
-                  IRNGT (IWRKD+3) = IRNG2
-!
-!   3 1 2
-!
-               Else
-                  IRNG1 = IRNGT (IWRKD+1)
-                  IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
-                  IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
-                  IRNGT (IWRKD+2) = IRNG1
-               End If
-               Exit
-            End If
-!
-!   1 2 3 4
-!
-            If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
-!
-!   1 3 x x
-!
-            If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
-               IRNG2 = IRNGT (IWRKD+2)
-               IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
-               If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
-!   1 3 2 4
-                  IRNGT (IWRKD+3) = IRNG2
-               Else
-!   1 3 4 2
-                  IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
-                  IRNGT (IWRKD+4) = IRNG2
-               End If
-!
-!   3 x x x
-!
-            Else
-               IRNG1 = IRNGT (IWRKD+1)
-               IRNG2 = IRNGT (IWRKD+2)
-               IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
-               If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
-                  IRNGT (IWRKD+2) = IRNG1
-                  If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
-!   3 1 2 4
-                     IRNGT (IWRKD+3) = IRNG2
-                  Else
-!   3 1 4 2
-                     IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
-                     IRNGT (IWRKD+4) = IRNG2
-                  End If
-               Else
-!   3 4 1 2
-                  IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
-                  IRNGT (IWRKD+3) = IRNG1
-                  IRNGT (IWRKD+4) = IRNG2
-               End If
-            End If
-         End Do
-!
-!  The Cs become As and Bs
-!
-         LMTNA = 4
-         Exit
-      End Do
-!
-!  Iteration loop. Each time, the length of the ordered subsets
-!  is doubled.
-!
-      Do
-         If (LMTNA >= NVAL) Exit
-         IWRKF = 0
-         LMTNC = 2 * LMTNC
-!
-!   Loop on merges of A and B into C
-!
-         Do
-            IWRK = IWRKF
-            IWRKD = IWRKF + 1
-            JINDA = IWRKF + LMTNA
-            IWRKF = IWRKF + LMTNC
-            If (IWRKF >= NVAL) Then
-               If (JINDA >= NVAL) Exit
-               IWRKF = NVAL
-            End If
-            IINDA = 1
-            IINDB = JINDA + 1
-!
-!   Shortcut for the case when the max of A is smaller
-!   than the min of B. This line may be activated when the
-!   initial set is already close to sorted.
-!
-!          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
-!
-!  One steps in the C subset, that we build in the final rank array
-!
-!  Make a copy of the rank array for the merge iteration
-!
-            JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
-!
-            XVALA = XDONT (JWRKT(IINDA))
-            XVALB = XDONT (IRNGT(IINDB))
-!
-            Do
-               IWRK = IWRK + 1
-!
-!  We still have unprocessed values in both A and B
-!
-               If (XVALA > XVALB) Then
-                  IRNGT (IWRK) = IRNGT (IINDB)
-                  IINDB = IINDB + 1
-                  If (IINDB > IWRKF) Then
-!  Only A still with unprocessed values
-                     IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
-                     Exit
-                  End If
-                  XVALB = XDONT (IRNGT(IINDB))
-               Else
-                  IRNGT (IWRK) = JWRKT (IINDA)
-                  IINDA = IINDA + 1
-                  If (IINDA > LMTNA) Exit! Only B still with unprocessed values
-                  XVALA = XDONT (JWRKT(IINDA))
-               End If
-!
-            End Do
-         End Do
-!
-!  The Cs become As and Bs
-!
-         LMTNA = 2 * LMTNA
-      End Do
-!
-      Return
-!
-End Subroutine I_mrgrnk
-end module m_mrgrnk
Index: LMDZ5/trunk/libf/phylmd/cosp/phys_cosp.F90.prev
===================================================================
--- LMDZ5/trunk/libf/phylmd/cosp/phys_cosp.F90.prev	(revision 2432)
+++ 	(revision )
@@ -1,456 +1,0 @@
-! Simulateur COSP : Cfmip Observation Simulator Package
-! ISCCP, Radar (QuickBeam), Lidar et Parasol (ACTSIM), MISR, RTTOVS
-!Idelkadi Abderrahmane Aout-Septembre 2009
-
-  subroutine phys_cosp( itap,dtime,freq_cosp,ecrit_mth,ecrit_day,ecrit_hf, &
-                        overlaplmdz,Nptslmdz,Nlevlmdz,lon,lat, presnivs, &
-                        ref_liq,ref_ice,fracTerLic,u_wind,v_wind,phi,ph,p,skt,t, &
-                        sh,rh,tca,cca,mr_lsliq,mr_lsice,fl_lsrainI,fl_lssnowI, &
-                        fl_ccrainI,fl_ccsnowI,mr_ozone,dtau_s,dem_s)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!! Inputs :
-! itap,                                 !Increment de la physiq
-! dtime,                                !Pas de temps physiq
-! overlap,                              !Overlap type in SCOPS
-! Npoints,                              !Nb de points de la grille physiq
-! Nlevels,                              !Nb de niveaux verticaux
-! Ncolumns,                             !Number of subcolumns
-! lon,lat,                              !Longitudes et latitudes de la grille LMDZ
-! ref_liq,ref_ice,                      !Rayons effectifs des particules liq et ice (en microm)
-! fracTerLic,                               !Fraction terre a convertir en masque
-! u_wind,v_wind,                        !Vents a 10m ???
-! phi,                                  !Geopotentiel
-! ph,                                   !pression pour chaque inter-couche
-! p,                                    !Pression aux milieux des couches
-! skt,t,                                !Temp au sol et temp 3D
-! sh,                                   !Humidite specifique
-! rh,                                   !Humidite relatif
-! tca,                                  !Fraction nuageuse
-! cca                                   !Fraction nuageuse convective
-! mr_lsliq,                             !Liq Cloud water content
-! mr_lsice,                             !Ice Cloud water content
-! mr_ccliq,                             !Convective Cloud Liquid water content  
-! mr_ccice,                             !Cloud ice water content
-! fl_lsrain,                            !Large scale precipitation lic
-! fl_lssnow,                            !Large scale precipitation ice
-! fl_ccrain,                            !Convective precipitation lic
-! fl_ccsnow,                            !Convective precipitation ice
-! mr_ozone,                             !Concentration ozone (Kg/Kg)
-! dem_s                                 !Cloud optical emissivity
-! dtau_s               			!Cloud optical thickness
-! emsfc_lw = 1.        			!Surface emissivity dans radlwsw.F90
-
-!!! Outputs :
-! calipso2D,                            !Lidar Low/heigh/Mean/Total-level Cloud Fraction
-! calipso3D,                            !Lidar Cloud Fraction (532 nm)
-! cfadlidar,                            !Lidar Scattering Ratio CFAD (532 nm)
-! parasolrefl,                          !PARASOL-like mono-directional reflectance
-! atb,                                  !Lidar Attenuated Total Backscatter (532 nm)
-! betamol,                              !Lidar Molecular Backscatter (532 nm)
-! cfaddbze,                             !Radar Reflectivity Factor CFAD (94 GHz)
-! clcalipso2,                           !Cloud frequency of occurrence as seen by CALIPSO but not CloudSat
-! dbze,                                 !Efective_reflectivity_factor
-! cltlidarradar,                        !Lidar and Radar Total Cloud Fraction
-! clMISR,                               !Cloud Fraction as Calculated by the MISR Simulator
-! clisccp2,                             !Cloud Fraction as Calculated by the ISCCP Simulator
-! boxtauisccp,                          !Optical Depth in Each Column as Calculated by the ISCCP Simulator
-! boxptopisccp,                         !Cloud Top Pressure in Each Column as Calculated by the ISCCP Simulator
-! tclisccp,                             !Total Cloud Fraction as Calculated by the ISCCP Simulator
-! ctpisccp,                             !Mean Cloud Top Pressure as Calculated by the ISCCP Simulator
-! tauisccp,                             !Mean Optical Depth as Calculated by the ISCCP Simulator
-! albisccp,                             !Mean Cloud Albedo as Calculated by the ISCCP Simulator
-! meantbisccp,                          !Mean all-sky 10.5 micron brightness temperature as calculated by the ISCCP Simulator
-! meantbclrisccp                        !Mean clear-sky 10.5 micron brightness temperature as calculated by the ISCCP Simulator
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-  USE MOD_COSP_CONSTANTS
-  USE MOD_COSP_TYPES
-  USE MOD_COSP
-  USE mod_phys_lmdz_para
-  use ioipsl
-  use iophy
- 
-  IMPLICIT NONE
-
-  ! Local variables
-  character(len=64)  :: cosp_input_nl='cosp_input_nl.txt'
-  character(len=64)  :: cosp_output_nl='cosp_output_nl.txt'
-  character(len=512), save :: finput ! Input file name
-  character(len=512), save :: cmor_nl
-  integer, save :: isccp_topheight,isccp_topheight_direction,overlap
-  integer,save  :: Ncolumns     ! Number of subcolumns in SCOPS
-  integer,parameter :: Ncollmdz=20
-  integer, save :: Npoints      ! Number of gridpoints
-  integer, save :: Nlevels      ! Number of levels
-  Integer :: Nptslmdz,Nlevlmdz ! Nb de points issus de physiq.F
-  integer, save :: Nlr          ! Number of levels in statistical outputs
-  integer, save :: Npoints_it   ! Max number of gridpoints to be processed in one iteration
-  integer :: i
-  type(cosp_config),save :: cfg   ! Configuration options
-  type(cosp_gridbox) :: gbx ! Gridbox information. Input for COSP
-  type(cosp_subgrid) :: sgx     ! Subgrid outputs
-  type(cosp_sgradar) :: sgradar ! Output from radar simulator
-  type(cosp_sglidar) :: sglidar ! Output from lidar simulator
-  type(cosp_isccp)   :: isccp   ! Output from ISCCP simulator
-  type(cosp_misr)    :: misr    ! Output from MISR simulator
-  type(cosp_vgrid)   :: vgrid   ! Information on vertical grid of stats
-  type(cosp_radarstats) :: stradar ! Summary statistics from radar simulator
-  type(cosp_lidarstats) :: stlidar ! Summary statistics from lidar simulator
-
-  integer :: t0,t1,count_rate,count_max
-  integer :: Nlon,Nlat,geomode
-  real,save :: radar_freq,k2,ZenAng,co2,ch4,n2o,co,emsfc_lw
-  integer,dimension(RTTOV_MAX_CHANNELS),save :: Channels
-  real,dimension(RTTOV_MAX_CHANNELS),save :: Surfem
-  integer, save :: surface_radar,use_mie_tables,use_gas_abs,do_ray,melt_lay
-  integer, save :: Nprmts_max_hydro,Naero,Nprmts_max_aero,lidar_ice_type
-  integer, save :: platform,satellite,Instrument,Nchannels
-  logical, save :: use_vgrid,csat_vgrid,use_precipitation_fluxes,use_reff
-
-! Declaration necessaires pour les sorties IOIPSL
-  integer :: ii,idayref
-  real    :: zjulian,zstoday,zstomth,zstohf,zout,ecrit_day,ecrit_hf,ecrit_mth
-  integer :: nhori,nvert,nvertp,nvertisccp,nvertm,nvertcol
-  integer, save :: nid_day_cosp,nid_mth_cosp,nid_hf_cosp
-  logical, save :: debut_cosp=.true.
-  integer :: itau_wcosp
-  character(len=10),dimension(Ncollmdz) :: chcol=(/'c01','c02','c03','c04','c05','c06','c07','c08','c09','c10', &
-                                                   'c11','c12','c13','c14','c15','c16','c17','c18','c19','c20'/)
-  real,dimension(Ncollmdz) :: column_ax
-  integer, save :: Nlevout
-
-  include "dimensions.h"
-  include "temps.h"  
-  
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Input variables from LMDZ-GCM
-  integer                         :: overlaplmdz   !  overlap type: 1=max, 2=rand, 3=max/rand ! cosp input (output lmdz)
-!  real,dimension(Npoints,Nlevels) :: height,phi,p,ph,T,sh,rh,tca,cca,mr_lsliq,mr_lsice,mr_ccliq,mr_ccice, &
-  real,dimension(Nptslmdz,Nlevlmdz) :: height,phi,p,ph,T,sh,rh,tca,cca,mr_lsliq,mr_lsice,mr_ccliq,mr_ccice, & 
-                                     fl_lsrain,fl_lssnow,fl_ccrain,fl_ccsnow,fl_lsgrpl, &
-                                     zlev,mr_ozone,radliq,radice,dtau_s,dem_s
-  real,dimension(Nptslmdz,Nlevlmdz) ::  fl_lsrainI,fl_lssnowI,fl_ccrainI,fl_ccsnowI
-  real,dimension(Nptslmdz)        :: lon,lat,skt,fracTerLic,u_wind,v_wind
-  real,dimension(Nlevlmdz)        :: presnivs
-  real                            :: ref_liq,ref_ice 
-  integer                         :: itap,k,ip
-  real                            :: dtime,freq_cosp
-
-!
-   namelist/COSP_INPUT/cmor_nl,overlap,isccp_topheight,isccp_topheight_direction, &
-              npoints,npoints_it,ncolumns,nlevels,use_vgrid,nlr,csat_vgrid,finput, &
-              radar_freq,surface_radar,use_mie_tables, &
-              use_gas_abs,do_ray,melt_lay,k2,Nprmts_max_hydro,Naero,Nprmts_max_aero, &
-              lidar_ice_type,use_precipitation_fluxes,use_reff, &
-              platform,satellite,Instrument,Nchannels, &
-              Channels,Surfem,ZenAng,co2,ch4,n2o,co
-
-!---------------- End of declaration of variables --------------
-   
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-! Read namelist with COSP inputs
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
- if (debut_cosp) then
-! Lecture du namelist input 
-  open(10,file=cosp_input_nl,status='old')
-  read(10,nml=cosp_input)
-  close(10)
-! Clefs Outputs 
-  call read_cosp_output_nl(cosp_output_nl,cfg)
-
-    if ( (Ncollmdz.ne.Ncolumns).or.(Nptslmdz.ne.Npoints).or.(Nlevlmdz.ne.Nlevels) ) then
-       print*,'Nb points Horiz, Vert, Sub-col passes par physiq.F = ', &
-               Nptslmdz, Nlevlmdz, Ncollmdz
-       print*,'Nb points Horiz, Vert, Sub-col lus dans namelist = ', &
-               Npoints, Nlevels, Ncolumns
-       print*,'Nb points Horiz, Vert, Sub-col passes par physiq.F est different de celui lu par namelist '
-       call abort
-    endif
-
-    if (overlaplmdz.ne.overlap) then
-       print*,'Attention overlaplmdz different de overlap lu dans namelist '
-    endif
-   print*,'Fin lecture Namelists, debut_cosp =',debut_cosp
-
-  endif ! debut_cosp
-
-  print*,'Debut phys_cosp itap,dtime,freq_cosp,ecrit_mth,ecrit_day,ecrit_hf ', &
-          itap,dtime,freq_cosp,ecrit_mth,ecrit_day,ecrit_hf
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-! Allocate local arrays
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-!        call system_clock(t0,count_rate,count_max) !!! Only for testing purposes
-        
-        
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-! Allocate memory for gridbox type
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        print *, 'Allocating memory for gridbox type...'
-        print*,'Npoints,Nlevels,Ncolumns,N_HYDRO,Nprmts_max_hydro ',Npoints,Nlevels,Ncolumns,N_HYDRO,Nprmts_max_hydro 
-  print*,' Cles sorties cosp :'
-  print*,' Lradar_sim,Llidar_sim,Lisccp_sim,Lmisr_sim,Lrttov_sim', &
-          cfg%Lradar_sim,cfg%Llidar_sim,cfg%Lisccp_sim,cfg%Lmisr_sim,cfg%Lrttov_sim
-        call construct_cosp_gridbox(float(itap),radar_freq,surface_radar,use_mie_tables,use_gas_abs,do_ray,melt_lay,k2, &
-                                    Npoints,Nlevels,Ncolumns,N_HYDRO,Nprmts_max_hydro,Naero,Nprmts_max_aero,Npoints_it, &
-                                    lidar_ice_type,isccp_topheight,isccp_topheight_direction,overlap,emsfc_lw, &
-                                    use_precipitation_fluxes,use_reff, &
-                                    Platform,Satellite,Instrument,Nchannels,ZenAng, &
-                                    channels(1:Nchannels),surfem(1:Nchannels),co2,ch4,n2o,co,gbx)
-        
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-! Here code to populate input structure
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-        print *, 'Populating input structure...'
-        gbx%longitude = lon
-        gbx%latitude = lat
-
-        gbx%p = p !
-        gbx%ph = ph
-        gbx%zlev_half = phi/9.81
-
-       do k = 1, Nlevels-1
-       do ip = 1, Npoints
-        zlev(ip,k) = phi(ip,k)/9.81 + (phi(ip,k+1)-phi(ip,k))/9.81 * (ph(ip,k)-ph(ip,k+1))/p(ip,k)
-       enddo
-       enddo
-       do ip = 1, Npoints
-        zlev(ip,Nlevels) = zlev(ip,Nlevels-1)+ 2.*(phi(ip,Nlevels)/9.81-zlev(ip,Nlevels-1))
-       END DO
-        gbx%zlev = zlev
-
-        gbx%T = T
-        gbx%q = rh*100.
-        gbx%sh = sh
-        gbx%cca = cca !convective_cloud_amount (1)
-        gbx%tca = tca ! total_cloud_amount (1)
-        gbx%psfc = ph(:,1) !pression de surface
-        gbx%skt  = skt !Skin temperature (K)
-
-        do ip = 1, Npoints
-          if (fracTerLic(ip).ge.0.5) then
-             gbx%land(ip) = 1.
-          else
-             gbx%land(ip) = 0.
-          endif
-        enddo
-        gbx%mr_ozone  = mr_ozone !mass_fraction_of_ozone_in_air (kg/kg)
-! A voir l equivalent LMDZ (u10m et v10m)
-        gbx%u_wind  = u_wind !eastward_wind (m s-1)
-        gbx%v_wind  = v_wind !northward_wind
-! Attention
-        gbx%sunlit  = 1
-
-! A voir l equivalent LMDZ
-  mr_ccliq = 0.0
-  mr_ccice = 0.0
-        gbx%mr_hydro(:,:,I_LSCLIQ) = mr_lsliq !mixing_ratio_large_scale_cloud_liquid (kg/kg)
-        gbx%mr_hydro(:,:,I_LSCICE) = mr_lsice !mixing_ratio_large_scale_cloud_ic
-        gbx%mr_hydro(:,:,I_CVCLIQ) = mr_ccliq !mixing_ratio_convective_cloud_liquid
-        gbx%mr_hydro(:,:,I_CVCICE) = mr_ccice !mixing_ratio_convective_cloud_ice
-! A revoir
-        fl_lsrain = fl_lsrainI + fl_ccrainI
-        fl_lssnow = fl_lssnowI + fl_ccsnowI
-        gbx%rain_ls = fl_lsrain !flux_large_scale_cloud_rain (kg m^-2 s^-1)
-        gbx%snow_ls = fl_lssnow !flux_large_scale_cloud_snow
-!  A voir l equivalent LMDZ
-        fl_lsgrpl=0.
-        fl_ccsnow = 0.
-        fl_ccrain = 0.
-        gbx%grpl_ls = fl_lsgrpl  !flux_large_scale_cloud_graupel
-        gbx%rain_cv = fl_ccrain  !flux_convective_cloud_rain
-        gbx%snow_cv = fl_ccsnow  !flux_convective_cloud_snow
-
-!Attention Teste
-       do k = 1, Nlevels
-        do ip = 1, Npoints
-!     liquid particles :
-         radliq(ip,k) = 12.0e-06
-         if (k.le.3) radliq(ip,k) = 11.0e-06
-
-!    ice particles :
-        if ( (t(ip,k)-273.15).gt.-81.4 ) then
-          radice(ip,k) = (0.71*(t(ip,k)-273.15)+61.29)*1e-6
-        else
-          radice(ip,k) = 3.5*1e-6
-        endif
-       END DO
-      END DO
-      gbx%Reff(:,:,I_LSCLIQ) = radliq
-      gbx%Reff(:,:,I_LSCICE) = radice
-      gbx%Reff(:,:,I_CVCLIQ) = radliq
-      gbx%Reff(:,:,I_CVCICE) = radice
-
-        ! ISCCP simulator
-        gbx%dtau_s   = dtau_s
-        print*,'dtau_s(1,:)=',gbx%dtau_s(1,:)
-        gbx%dtau_c   = 0.
-        gbx%dem_s    = dem_s
-        print*,'dem_s(1,:)=',gbx%dem_s(1,:)
-        gbx%dem_c    = 0.
-
-! Surafce emissivity
-       emsfc_lw = 1.
-               
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        ! Define new vertical grid
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        print *, 'Defining new vertical grid...'
-        call construct_cosp_vgrid(gbx,Nlr,use_vgrid,csat_vgrid,vgrid)
-
- if (debut_cosp) then
-! Creer le fichier de sorie, definir les variable de sortie
-  ! Axe verticale (Pa ou Km)
-     Nlevout = vgrid%Nlvgrid
-   
-        do ii=1,Ncolumns
-          column_ax(ii) = float(ii)
-        enddo
-
-     include "ini_histmthCOSP.h"
-     include "ini_histdayCOSP.h"
-     include "ini_histhfCOSP.h"
-
-   print*,'Fin Initialisation des sorties COSP, debut_cosp =',debut_cosp 
-   print*,'R_UNDEF=',R_UNDEF
-
-   debut_cosp=.false.
-  endif ! debut_cosp
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-       ! Allocate memory for other types
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        print *, 'Allocating memory for other types...'
-        call construct_cosp_subgrid(Npoints, Ncolumns, Nlevels, sgx)
-        call construct_cosp_sgradar(cfg,Npoints,Ncolumns,Nlevels,N_HYDRO,sgradar)
-        call construct_cosp_radarstats(cfg,Npoints,Ncolumns,vgrid%Nlvgrid,N_HYDRO,stradar)
-        call construct_cosp_sglidar(cfg,Npoints,Ncolumns,Nlevels,N_HYDRO,PARASOL_NREFL,sglidar)
-        call construct_cosp_lidarstats(cfg,Npoints,Ncolumns,vgrid%Nlvgrid,N_HYDRO,PARASOL_NREFL,stlidar)
-        call construct_cosp_isccp(cfg,Npoints,Ncolumns,Nlevels,isccp)
-        call construct_cosp_misr(cfg,Npoints,misr)
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        ! Call simulator
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        print *, 'Calling simulator...'
-        call cosp(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,stradar,stlidar)
-         print*,'stlidar%lidarcld(1,:)=',stlidar%lidarcld(1,:)
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        ! Write outputs to CMOR-compliant NetCDF
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-! A traiter le cas ou l on a des valeurs indefinies
-! Attention teste
-
-! if(1.eq.0)then
-
-
-   do k = 1,Nlevout
-!     do ip = 1,Npoints
-!     if(stlidar%lidarcld(ip,k).eq.R_UNDEF)then
-!      stlidar%lidarcld(ip,k)=0.
-!     endif
-!     enddo
-
-
-     do ii= 1,SR_BINS
-      do ip = 1,Npoints
-       if(stlidar%cfad_sr(ip,ii,k).eq.R_UNDEF)then
-        stlidar%cfad_sr(ip,ii,k)=0.
-       endif
-      enddo
-     enddo
-   enddo   
-   
-  do ip = 1,Npoints
-   do k = 1,Nlevlmdz 
-     if(sglidar%beta_mol(ip,k).eq.R_UNDEF)then
-      sglidar%beta_mol(ip,k)=0.
-     endif
-    
-     do ii= 1,Ncolumns
-       if(sglidar%beta_tot(ip,ii,k).eq.R_UNDEF)then
-        sglidar%beta_tot(ip,ii,k)=0.
-       endif  
-     enddo
-
-    enddo    !k = 1,Nlevlmdz
-   enddo     !ip = 1,Npoints
-
-   do k = 1,LIDAR_NCAT
-    do ip = 1,Npoints
-     if(stlidar%cldlayer(ip,k).eq.R_UNDEF)then
-      stlidar%cldlayer(ip,k)=0.
-     endif
-    enddo
-   enddo
-
-! endif 
-
-   do ip = 1,Npoints
-    if(isccp%totalcldarea(ip).eq.-1.E+30)then
-      isccp%totalcldarea(ip)=0.
-    endif
-    if(isccp%meanptop(ip).eq.-1.E+30)then
-      isccp%meanptop(ip)=0.
-    endif
-    if(isccp%meantaucld(ip).eq.-1.E+30)then
-      isccp%meantaucld(ip)=0.
-    endif
-    if(isccp%meanalbedocld(ip).eq.-1.E+30)then
-      isccp%meanalbedocld(ip)=0.
-    endif
-    if(isccp%meantb(ip).eq.-1.E+30)then
-      isccp%meantb(ip)=0.
-    endif
-    if(isccp%meantbclr(ip).eq.-1.E+30)then
-      isccp%meantbclr(ip)=0.
-    endif
-
-    do k=1,7
-     do ii=1,7
-     if(isccp%fq_isccp(ip,ii,k).eq.-1.E+30)then
-      isccp%fq_isccp(ip,ii,k)=0.
-     endif
-     enddo
-    enddo
-
-    do ii=1,Ncolumns
-     if(isccp%boxtau(ip,ii).eq.-1.E+30)then
-       isccp%boxtau(ip,ii)=0.
-     endif
-    enddo
-
-    do ii=1,Ncolumns
-     if(isccp%boxptop(ip,ii).eq.-1.E+30)then
-       isccp%boxptop(ip,ii)=0.
-     endif
-    enddo
-   enddo
-
-  include "write_histmthCOSP.h"
-  include "write_histdayCOSP.h"
-  include "write_histhfCOSP.h"
-
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        ! Deallocate memory in derived types
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        print *, 'Deallocating memory...'
-        call free_cosp_gridbox(gbx)
-        call free_cosp_subgrid(sgx)
-        call free_cosp_sgradar(sgradar)
-        call free_cosp_radarstats(stradar)
-        call free_cosp_sglidar(sglidar)
-        call free_cosp_lidarstats(stlidar)
-        call free_cosp_isccp(isccp)
-        call free_cosp_misr(misr)
-        call free_cosp_vgrid(vgrid)  
-  
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  ! Time in s. Only for testing purposes
-!  call system_clock(t1,count_rate,count_max)
-!  print *,(t1-t0)*1.0/count_rate
-    
-end subroutine phys_cosp
Index: LMDZ5/trunk/libf/phylmd/cosp/test_modis_simulator.F90.s
===================================================================
--- LMDZ5/trunk/libf/phylmd/cosp/test_modis_simulator.F90.s	(revision 2432)
+++ 	(revision )
@@ -1,154 +1,0 @@
-! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $
-! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/MODIS_simulator/test_modis_simulator.F90 $
-program test_modis_simulator
-  use mod_modis_sim
-  implicit none
-  ! 
-  ! Tests cases for MODIS L2 (pixel) and L3 (grid-scale) simulators
-  !
-  
-  !
-  ! Inputs
-  !
-  integer, parameter :: nLayers = 6, nSubCols = 10
-  real, dimension(1, nLayers)    :: temp, pressureLayers
-  real, dimension(1, nLayers+1)  :: pressureLevels
-  real, dimension(1, nSubCols, nLayers) :: &
-                                    opticalThickness, cloudWater, cloudIce, waterSize, iceSize, &
-                                    liquid_opticalThickness, ice_opticalThickness
-  !
-  ! Pixel-scale retreivals
-  ! 
-  integer, dimension(1, nSubCols) :: retrievedPhase
-  real,    dimension(1, nSubCols) :: isccpTau, isccpCloudTopPressure, &
-                                  retrievedCloudTopPressure, retrievedTau, retrievedSize
-  integer :: i, k
-  character(len = 6) :: phase
-  
-  !
-  ! Grid mean properties 
-  !
-  real, dimension(1) :: &
-    Cloud_Fraction_Total_Mean,       Cloud_Fraction_Water_Mean,       Cloud_Fraction_Ice_Mean,       &
-    Optical_Thickness_Total_Mean,    Optical_Thickness_Water_Mean,    Optical_Thickness_Ice_Mean,    &
-    Optical_Thickness_Total_LogMean, Optical_Thickness_Water_LogMean, Optical_Thickness_Ice_LogMean, &
-                                     Cloud_Particle_Size_Water_Mean,  Cloud_Particle_Size_Ice_Mean,  &
-    Cloud_Top_Pressure_Total_Mean,                                                                   &
-                                     Liquid_Water_Path_Mean,          Ice_Water_Path_Mean
-  real, dimension(1, numTauHistogramBins, numPressureHistogramBins) :: Optical_Thickness_vs_Cloud_Top_Pressure
-
-  ! ---------------------------------------------------------------------------------------------
-  pressureLevels(1, :) = (/   0., 200., 400., 600., 800., 1000., 1200. /) 
-  pressureLayers(1, :) = (pressureLevels(1, 2:) + pressureLevels(1, 1:nLayers)) / 2
-  temp(1, :) =           (/ 200., 220., 240., 260., 280., 300.0 /) 
-  isccpCloudTopPressure(:, :) = 999. 
-  
-
-  !
-  ! Tests:  these 9 from Steve Platnick
-  !
-  ice_opticalThickness = 0.; liquid_opticalThickness = 0.
-  iceSize              = 0.; waterSize = 0. 
-  
-  ice_opticalThickness   (1, 1, 3) = 0.2
-  iceSize                (1, 1, 3) = 30
-  
-  ice_opticalThickness   (1, 2, 3) = 1.
-  iceSize                (1, 2, 3) = 30
-  
-  liquid_opticalThickness(1, 3, 5) = 0.2
-  waterSize              (1, 3, 5) = 10 
-
-  liquid_opticalThickness(1, 4, 5) = 1.
-  waterSize              (1, 4, 5) = 10 
-  
-  ice_opticalThickness   (1, 5, 3) = 0.2
-  iceSize                (1, 5, 3) = 30
-  liquid_opticalThickness(1, 5, 5) = 5.
-  waterSize              (1, 5, 5) = 10 
-
-  ice_opticalThickness   (1, 6, 3) = 1. 
-  iceSize                (1, 6, 3) = 30
-  liquid_opticalThickness(1, 6, 5) = 5.
-  waterSize              (1, 6, 5) = 10
-  
-  ice_opticalThickness   (1, 7, 3) = 10. 
-  iceSize                (1, 7, 3) = 30
-  liquid_opticalThickness(1, 7, 5) = 5.
-  waterSize              (1, 7, 5) = 10
-  
-  ice_opticalThickness   (1, 8, 2) = 1. 
-  iceSize                (1, 8, 2) = 15.
-  ice_opticalThickness   (1, 8, 3) = 1. 
-  iceSize                (1, 8, 3) = 30.
-  ice_opticalThickness   (1, 8, 4) = 1. 
-  iceSize                (1, 8, 4) = 60.
-  
-  liquid_opticalThickness(1, 9, 4) = 1. 
-  waterSize              (1, 9, 4) = 15
-  liquid_opticalThickness(1, 9, 5) = 3. 
-  waterSize              (1, 9, 5) = 10.
-  liquid_opticalThickness(1, 9, 6) = 1. 
-  waterSize              (1, 9, 6) = 50.
-  
-  call modis_l2_simulator(temp(1,:), pressureLayers(1,:), pressureLevels(1,:),                       &
-                         liquid_opticalThickness(1,:, :), ice_opticalThickness(1,:, :), &
-                         waterSize(1,:, :), iceSize(1,:, :), &
-                         isccpTau(1, :), isccpCloudTopPressure(1, :),                            &
-                         retrievedPhase(1, :), retrievedCloudTopPressure(1, :), retrievedTau(1, :), retrievedSize(1, :))
-  do i = 1, nSubcols
-    print * 
-    print *, "Profile ", i
-    if(any(ice_opticalThickness(1, i, :) > 0.) .or. any(liquid_opticalThickness(1, i, :) > 0.)) & 
-      print *, "center p  ice: tau   size    water: tau    size" 
-    do k = 1, nLayers
-      if(ice_opticalThickness(1, i, k) > 0. .or. liquid_opticalThickness(1, i, k) > 0.) &
-        write(*, '(3x, f5.0, 7x, 2(f4.1, 3x, f4.1, 10x))'), &
-          pressureLayers(1, k), ice_opticalThickness(1, i, k), iceSize(1, i, k), liquid_opticalThickness(1, i, k), waterSize(1, i, k)
-    end do
-    select case(retrievedPhase(1, i)) 
-      case(phaseIsNone) 
-        phase = "None" 
-      case(phaseIsIce) 
-        phase = "Ice" 
-      case(phaseIsLiquid)
-        phase = "Liquid" 
-      case(phaseIsUndetermined)
-        phase = "Undet."
-      case default
-        phase = "???" 
-    end select
-    print *, "Retrievals" 
-    write(*, '("  Phase: ", a6, 2x, "CTP: ", f5.0, ", tau = ", f4.1, ", size = ", f4.1)') &
-          phase, retrievedCloudTopPressure(1, i), retrievedTau(1, i), retrievedSize(1, i)
-  end do 
-  
-  call modis_L3_simulator(retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize,            &
-       Cloud_Fraction_Total_Mean,       Cloud_Fraction_Water_Mean,       Cloud_Fraction_Ice_Mean,       &
-       Optical_Thickness_Total_Mean,    Optical_Thickness_Water_Mean,    Optical_Thickness_Ice_Mean,    &
-       Optical_Thickness_Total_LogMean, Optical_Thickness_Water_LogMean, Optical_Thickness_Ice_LogMean, &
-                                        Cloud_Particle_Size_Water_Mean,  Cloud_Particle_Size_Ice_Mean,  &
-       Cloud_Top_Pressure_Total_Mean,                                                                   &
-                                        Liquid_Water_Path_Mean,          Ice_Water_Path_Mean,           &    
-       Optical_Thickness_vs_Cloud_Top_Pressure)
-
-  print *; print *; print *, "Grid means" 
-  do i = 1, nSubcols
-    write(*, '("  Phase: ", i2, 2x, "CTP: ", f5.0, ", tau = ", f4.1, ", size = ", f4.1)') &
-          retrievedPhase(1, i), retrievedCloudTopPressure(1, i), retrievedTau(1, i), retrievedSize(1, i)
-  end do
-  print *, "         Total    Liquid     Ice" 
-  write(*, '("Fraction: ", 3(f4.3, 6x))') Cloud_Fraction_Total_Mean,       Cloud_Fraction_Water_Mean,       Cloud_Fraction_Ice_Mean
-  write(*, '("Tau:      ", 3(f4.1, 6x))')  Optical_Thickness_Total_Mean,    Optical_Thickness_Water_Mean,    Optical_Thickness_Ice_Mean
-  write(*, '("TauLog:   ", 3(f4.1, 6x))')  Optical_Thickness_Total_LogMean, Optical_Thickness_Water_LogMean, Optical_Thickness_Ice_LogMean
-  write(*, '("Size:     ", 10x, 2(f4.1, 6x))')  Cloud_Particle_Size_Water_Mean,  Cloud_Particle_Size_Ice_Mean
-  write(*, '("WaterPath:",  9x, 2(f5.1, 5x))')  Liquid_Water_Path_Mean,          Ice_Water_Path_Mean
-  write(*, '("CTP:      ", f5.1)') Cloud_Top_Pressure_Total_Mean
-  print *; print *, "Histogram" 
-  write(*, '(8x, 7(f4.1, 3x))')  tauHistogramBoundaries(:)
-  do k = 1, numPressureHistogramBins
-     write(*, '(f5.0, 3x, 7(f4.2, 3x))') pressureHistogramBoundaries(k), Optical_Thickness_vs_Cloud_Top_Pressure(1, :, k)
-  end do 
-  print *, "Total fraction from histogram:", sum(Optical_Thickness_vs_Cloud_Top_Pressure(1, :, :))
-end program   test_modis_simulator
-
