Index: LMDZ5/trunk/libf/phylmd/readaerosol_interp.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/readaerosol_interp.F90 (revision 2838) +++ 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 !****************************************************************************************