Index: /LMDZ5/trunk/libf/phylmd/readaerosol_interp.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/readaerosol_interp.F90	(revision 2839)
+++ /LMDZ5/trunk/libf/phylmd/readaerosol_interp.F90	(revision 2840)
@@ -162,5 +162,5 @@
      NULLIFY(pt_ap)
      NULLIFY(pt_b)
-  END IF
+  ENDIF
 
 !****************************************************************************************
@@ -187,5 +187,5 @@
            filename='aerosols'
            type='annuel'
-        END IF
+        ENDIF
      ELSE  IF (aer_type == 'mix2') THEN
         ! Special case using a mix of decenal sulfate file and natrual aerosols
@@ -196,5 +196,5 @@
            filename='aerosols'
            type='preind'
-        END IF
+        ENDIF
      ELSE  IF (aer_type == 'mix3') THEN
         ! Special case using a mix of annual sulfate file and natrual aerosols
@@ -205,8 +205,8 @@
            filename='aerosols' 
            type='preind'
-        END IF
+        ENDIF
      ELSE
         CALL abort_physic('readaerosol_interp', 'this aer_type not supported',1)
-     END IF
+     ENDIF
 
      CALL readaerosol(name_aero(id_aero), type, filename, iyr, klev_src, pt_ap, pt_b, pt_tmp, &
@@ -215,5 +215,5 @@
         ALLOCATE(var_year(klon, klev_src, 12, naero_spc), stat=ierr)
         IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 5',1)
-     END IF
+     ENDIF
      var_year(:,:,:,id_aero) = pt_tmp(:,:,:)
 
@@ -229,10 +229,10 @@
         WRITE(lunout,*) 'Aerosol : ', name_aero(id_aero)
         CALL abort_physic('readaerosol_interp','Differnt vertical axes in aerosol forcing files',1)
-     END IF
+     ENDIF
 
      IF (.NOT. ALLOCATED(pi_var_year)) THEN
         ALLOCATE(pi_var_year(klon, klev_src, 12, naero_spc), stat=ierr)
         IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 6',1)
-     END IF
+     ENDIF
      pi_var_year(:,:,:,id_aero) = pt_tmp(:,:,:)
     
@@ -244,5 +244,5 @@
         CALL writefield_phy('load_year_src',load_year(:,:,id_aero),1)
         CALL writefield_phy('pi_load_year_src',pi_load_year(:,:,id_aero),1)
-     END IF
+     ENDIF
 
      ! Pointer no more useful, deallocate. 
@@ -258,14 +258,14 @@
            WRITE(lunout,*) 'Warning! All forcing files for the same aerosol must have the same structure'
            CALL abort_physic('readaerosol_interp', 'The aerosol files have not the same format',1)
-        END IF
+        ENDIF
         
         IF (klev /= klev_src) THEN
            WRITE(lunout,*) 'Old format of aerosol file do not allowed vertical interpolation'
            CALL abort_physic('readaerosol_interp', 'Old aerosol file not possible',1)
-        END IF
+        ENDIF
 
      ELSE 
         vert_interp = .TRUE.
-     END IF
+     ENDIF
 
 !    Calendar initialisation
@@ -286,5 +286,5 @@
      endif
 
-  END IF  ! IF ( (first .OR. iday==0) .AND. lnewday ) THEN 
+  ENDIF  ! IF ( (first .OR. iday==0) .AND. lnewday ) THEN 
   
 !****************************************************************************************
@@ -309,5 +309,5 @@
              ! the month is january, thus the month before december
              im2=12
-          END IF
+          ENDIF
        ELSE
           ! the second half of the month
@@ -319,5 +319,5 @@
              im2=1
           ENDIF
-       END IF
+       ENDIF
      ELSE IF (nbr_tsteps == 14) then
        im = im + 1
@@ -332,5 +332,5 @@
           day1 = month_mid(im)
           day2 = month_mid(im2)
-       END IF
+       ENDIF
      ELSE
        CALL abort_physic('readaerosol_interp', 'number of months undefined',1)
@@ -358,6 +358,6 @@
                 pi_var_year(i,k,im2,id_aero) - (jDay-day2)/(day1-day2) * &
                 (pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero))
-        END DO
-     END DO
+        ENDDO
+     ENDDO
 
      ! Time interpolation for pressure at surface, still on vertical source grid
@@ -370,5 +370,5 @@
              pi_psurf_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * &
              (pi_psurf_year(i,im2,id_aero) - pi_psurf_year(i,im,id_aero))
-     END DO
+     ENDDO
 
      ! Time interpolation for the load, still on vertical source grid
@@ -381,5 +381,5 @@
              pi_load_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * &
              (pi_load_year(i,im2,id_aero) - pi_load_year(i,im,id_aero))
-     END DO
+     ENDDO
 
 !****************************************************************************************
@@ -397,6 +397,6 @@
            DO i = 1, klon
               pplay_src(i,k)= pt_ap(k) + pt_b(k)*psurf_day(i)
-           END DO
-        END DO
+           ENDDO
+        ENDDO
         
         IF (debug) THEN
@@ -406,5 +406,5 @@
            CALL writefield_phy('day_src',tmp1,klev_src)
            CALL writefield_phy('pi_day_src',tmp2,klev_src)
-        END IF
+        ENDIF
 
         ! b) vertical interpolation on pressure leveles
@@ -422,6 +422,6 @@
            DO i = 1, klon
               delp(i,k) = paprs(i,k) - paprs (i,k+1)
-           END DO
-        END DO
+           ENDDO
+        ENDDO
 
         ! Find the mass load in the actual pillar, on target grid
@@ -431,14 +431,14 @@
               zrho = pplay(i,k)/t_seri(i,k)/RD       ! [kg/m3]
               volm = var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
-              load_tgt(i) = load_tgt(i) + 1/RG * volm *delp(i,k)
-           END DO
-        END DO
+              load_tgt(i) = load_tgt(i) + volm *delp(i,k)/RG
+           ENDDO
+        ENDDO
         
         ! Adjust, uniform
         DO k = 1, klev
            DO i = 1, klon
-              var_day(i,k,id_aero) = var_day(i,k,id_aero)*load_src(i)/load_tgt(i) 
-           END DO
-        END DO
+              var_day(i,k,id_aero) = var_day(i,k,id_aero)*load_src(i)/max(1.e-30,load_tgt(i))
+           ENDDO
+        ENDDO
         
         IF (debug) THEN
@@ -448,7 +448,7 @@
                  zrho = pplay(i,k)/t_seri(i,k)/RD       ! [kg/m3]
                  volm = var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
-                 load_tgt_test(i) = load_tgt_test(i) + 1/RG * volm*delp(i,k)
-              END DO
-           END DO
+                 load_tgt_test(i) = load_tgt_test(i) + volm*delp(i,k)/RG
+              ENDDO
+           ENDDO
            
            CALL writefield_phy('day_tgt2',var_day(:,:,id_aero),klev)
@@ -456,5 +456,5 @@
            CALL writefield_phy('load_tgt_test',load_tgt_test(:),1)
            CALL writefield_phy('load_src',load_src(:),1)
-        END IF
+        ENDIF
 
         ! - Interpolate variable tmp2 (source grid) to pi_var_day (target grid)
@@ -464,11 +464,11 @@
            DO i = 1, klon
               pplay_src(i,k)= pt_ap(k) + pt_b(k)*pi_psurf_day(i)
-           END DO
-        END DO
+           ENDDO
+        ENDDO
 
         IF (debug) THEN
            CALL writefield_phy('pi_psurf_day_src',pi_psurf_day(:),1)
            CALL writefield_phy('pi_pplay_src',pplay_src(:,:),klev_src)
-        END IF
+        ENDIF
 
         ! b) vertical interpolation on pressure leveles
@@ -488,13 +488,13 @@
               zrho = pplay(i,k)/t_seri(i,k)/RD          ! [kg/m3]
               volm = pi_var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
-              load_tgt(i) = load_tgt(i) + 1/RG * volm * delp(i,k)
-           END DO
-        END DO
+              load_tgt(i) = load_tgt(i) + volm*delp(i,k)/RG
+           ENDDO
+        ENDDO
 
         DO k = 1, klev
            DO i = 1, klon
-              pi_var_day(i,k,id_aero) = pi_var_day(i,k,id_aero)*pi_load_src(i)/load_tgt(i)
-           END DO
-        END DO
+              pi_var_day(i,k,id_aero) = pi_var_day(i,k,id_aero)*pi_load_src(i)/max(1.e-30,load_tgt(i))
+           ENDDO
+        ENDDO
 
         IF (debug) THEN
@@ -504,12 +504,12 @@
                  zrho = pplay(i,k)/t_seri(i,k)/RD          ! [kg/m3]
                  volm = pi_var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
-                 load_tgt_test(i) = load_tgt_test(i) + 1/RG * volm * delp(i,k)
-              END DO
-           END DO
+                 load_tgt_test(i) = load_tgt_test(i) + volm*delp(i,k)/RG
+              ENDDO
+           ENDDO
            CALL writefield_phy('pi_day_tgt2',pi_var_day(:,:,id_aero),klev)
            CALL writefield_phy('pi_load_tgt',load_tgt(:),1)
            CALL writefield_phy('pi_load_tgt_test',load_tgt_test(:),1)
            CALL writefield_phy('pi_load_src',pi_load_src(:),1)
-        END IF
+        ENDIF
 
 
@@ -519,5 +519,5 @@
         pi_var_day(:,:,id_aero) = tmp2(:,:)
 
-     END IF ! vert_interp
+     ENDIF ! vert_interp
 
 
@@ -539,12 +539,12 @@
                          trim(name_aero(id_aero)),'(i,k,im)=',           &
                          var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero)
-                 END IF
+                 ENDIF
                  WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero)
                  WRITE(lunout,*) 'day1, day2, jDay = ', day1, day2, jDay 
                  CALL abort_physic('readaerosol_interp','Error in interpolation 1',1)
-              END IF
-           END DO 
-        END DO
-     END IF
+              ENDIF
+           ENDDO 
+        ENDDO
+     ENDIF
 
      IF (MINVAL(pi_var_day(:,:,id_aero)) < 0. ) THEN
@@ -558,14 +558,14 @@
                          trim(name_aero(id_aero)),'(i,k,im)=',           &
                          pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero)
-                 END IF
+                 ENDIF
                  
                  WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero)
                  CALL abort_physic('readaerosol_interp','Error in interpolation 2',1)
-              END IF
-           END DO
-        END DO
-     END IF
-
-  END IF ! lnewday
+              ENDIF
+           ENDDO
+        ENDDO
+     ENDIF
+
+  ENDIF ! lnewday
 
 !****************************************************************************************
