Changeset 2840


Ignore:
Timestamp:
Mar 30, 2017, 6:28:00 PM (8 years ago)
Author:
oboucher
Message:

Correcting the case when aerosol is zero (eg NO3 when missing from file)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/readaerosol_interp.F90

    r2523 r2840  
    162162     NULLIFY(pt_ap)
    163163     NULLIFY(pt_b)
    164   END IF
     164  ENDIF
    165165
    166166!****************************************************************************************
    … …  
    187187           filename='aerosols'
    188188           type='annuel'
    189         END IF
     189        ENDIF
    190190     ELSE  IF (aer_type == 'mix2') THEN
    191191        ! Special case using a mix of decenal sulfate file and natrual aerosols
    … …  
    196196           filename='aerosols'
    197197           type='preind'
    198         END IF
     198        ENDIF
    199199     ELSE  IF (aer_type == 'mix3') THEN
    200200        ! Special case using a mix of annual sulfate file and natrual aerosols
    … …  
    205205           filename='aerosols'
    206206           type='preind'
    207         END IF
     207        ENDIF
    208208     ELSE
    209209        CALL abort_physic('readaerosol_interp', 'this aer_type not supported',1)
    210      END IF
     210     ENDIF
    211211
    212212     CALL readaerosol(name_aero(id_aero), type, filename, iyr, klev_src, pt_ap, pt_b, pt_tmp, &
    … …  
    215215        ALLOCATE(var_year(klon, klev_src, 12, naero_spc), stat=ierr)
    216216        IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 5',1)
    217      END IF
     217     ENDIF
    218218     var_year(:,:,:,id_aero) = pt_tmp(:,:,:)
    219219
    … …  
    229229        WRITE(lunout,*) 'Aerosol : ', name_aero(id_aero)
    230230        CALL abort_physic('readaerosol_interp','Differnt vertical axes in aerosol forcing files',1)
    231      END IF
     231     ENDIF
    232232
    233233     IF (.NOT. ALLOCATED(pi_var_year)) THEN
    234234        ALLOCATE(pi_var_year(klon, klev_src, 12, naero_spc), stat=ierr)
    235235        IF (ierr /= 0) CALL abort_physic('readaerosol_interp', 'pb in allocation 6',1)
    236      END IF
     236     ENDIF
    237237     pi_var_year(:,:,:,id_aero) = pt_tmp(:,:,:)
    238238   
    … …  
    244244        CALL writefield_phy('load_year_src',load_year(:,:,id_aero),1)
    245245        CALL writefield_phy('pi_load_year_src',pi_load_year(:,:,id_aero),1)
    246      END IF
     246     ENDIF
    247247
    248248     ! Pointer no more useful, deallocate.
    … …  
    258258           WRITE(lunout,*) 'Warning! All forcing files for the same aerosol must have the same structure'
    259259           CALL abort_physic('readaerosol_interp', 'The aerosol files have not the same format',1)
    260         END IF
     260        ENDIF
    261261       
    262262        IF (klev /= klev_src) THEN
    263263           WRITE(lunout,*) 'Old format of aerosol file do not allowed vertical interpolation'
    264264           CALL abort_physic('readaerosol_interp', 'Old aerosol file not possible',1)
    265         END IF
     265        ENDIF
    266266
    267267     ELSE
    268268        vert_interp = .TRUE.
    269      END IF
     269     ENDIF
    270270
    271271!    Calendar initialisation
    … …  
    286286     endif
    287287
    288   END IF  ! IF ( (first .OR. iday==0) .AND. lnewday ) THEN
     288  ENDIF  ! IF ( (first .OR. iday==0) .AND. lnewday ) THEN
    289289 
    290290!****************************************************************************************
    … …  
    309309             ! the month is january, thus the month before december
    310310             im2=12
    311           END IF
     311          ENDIF
    312312       ELSE
    313313          ! the second half of the month
    … …  
    319319             im2=1
    320320          ENDIF
    321        END IF
     321       ENDIF
    322322     ELSE IF (nbr_tsteps == 14) then
    323323       im = im + 1
    … …  
    332332          day1 = month_mid(im)
    333333          day2 = month_mid(im2)
    334        END IF
     334       ENDIF
    335335     ELSE
    336336       CALL abort_physic('readaerosol_interp', 'number of months undefined',1)
    … …  
    358358                pi_var_year(i,k,im2,id_aero) - (jDay-day2)/(day1-day2) * &
    359359                (pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero))
    360         END DO
    361      END DO
     360        ENDDO
     361     ENDDO
    362362
    363363     ! Time interpolation for pressure at surface, still on vertical source grid
    … …  
    370370             pi_psurf_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * &
    371371             (pi_psurf_year(i,im2,id_aero) - pi_psurf_year(i,im,id_aero))
    372      END DO
     372     ENDDO
    373373
    374374     ! Time interpolation for the load, still on vertical source grid
    … …  
    381381             pi_load_year(i,im2,id_aero) - (jDay-day2)/(day1-day2) * &
    382382             (pi_load_year(i,im2,id_aero) - pi_load_year(i,im,id_aero))
    383      END DO
     383     ENDDO
    384384
    385385!****************************************************************************************
    … …  
    397397           DO i = 1, klon
    398398              pplay_src(i,k)= pt_ap(k) + pt_b(k)*psurf_day(i)
    399            END DO
    400         END DO
     399           ENDDO
     400        ENDDO
    401401       
    402402        IF (debug) THEN
    … …  
    406406           CALL writefield_phy('day_src',tmp1,klev_src)
    407407           CALL writefield_phy('pi_day_src',tmp2,klev_src)
    408         END IF
     408        ENDIF
    409409
    410410        ! b) vertical interpolation on pressure leveles
    … …  
    422422           DO i = 1, klon
    423423              delp(i,k) = paprs(i,k) - paprs (i,k+1)
    424            END DO
    425         END DO
     424           ENDDO
     425        ENDDO
    426426
    427427        ! Find the mass load in the actual pillar, on target grid
    … …  
    431431              zrho = pplay(i,k)/t_seri(i,k)/RD       ! [kg/m3]
    432432              volm = var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
    433               load_tgt(i) = load_tgt(i) + 1/RG * volm *delp(i,k)
    434            END DO
    435         END DO
     433              load_tgt(i) = load_tgt(i) + volm *delp(i,k)/RG
     434           ENDDO
     435        ENDDO
    436436       
    437437        ! Adjust, uniform
    438438        DO k = 1, klev
    439439           DO i = 1, klon
    440               var_day(i,k,id_aero) = var_day(i,k,id_aero)*load_src(i)/load_tgt(i)
    441            END DO
    442         END DO
     440              var_day(i,k,id_aero) = var_day(i,k,id_aero)*load_src(i)/max(1.e-30,load_tgt(i))
     441           ENDDO
     442        ENDDO
    443443       
    444444        IF (debug) THEN
    … …  
    448448                 zrho = pplay(i,k)/t_seri(i,k)/RD       ! [kg/m3]
    449449                 volm = var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
    450                  load_tgt_test(i) = load_tgt_test(i) + 1/RG * volm*delp(i,k)
    451               END DO
    452            END DO
     450                 load_tgt_test(i) = load_tgt_test(i) + volm*delp(i,k)/RG
     451              ENDDO
     452           ENDDO
    453453           
    454454           CALL writefield_phy('day_tgt2',var_day(:,:,id_aero),klev)
    … …  
    456456           CALL writefield_phy('load_tgt_test',load_tgt_test(:),1)
    457457           CALL writefield_phy('load_src',load_src(:),1)
    458         END IF
     458        ENDIF
    459459
    460460        ! - Interpolate variable tmp2 (source grid) to pi_var_day (target grid)
    … …  
    464464           DO i = 1, klon
    465465              pplay_src(i,k)= pt_ap(k) + pt_b(k)*pi_psurf_day(i)
    466            END DO
    467         END DO
     466           ENDDO
     467        ENDDO
    468468
    469469        IF (debug) THEN
    470470           CALL writefield_phy('pi_psurf_day_src',pi_psurf_day(:),1)
    471471           CALL writefield_phy('pi_pplay_src',pplay_src(:,:),klev_src)
    472         END IF
     472        ENDIF
    473473
    474474        ! b) vertical interpolation on pressure leveles
    … …  
    488488              zrho = pplay(i,k)/t_seri(i,k)/RD          ! [kg/m3]
    489489              volm = pi_var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
    490               load_tgt(i) = load_tgt(i) + 1/RG * volm * delp(i,k)
    491            END DO
    492         END DO
     490              load_tgt(i) = load_tgt(i) + volm*delp(i,k)/RG
     491           ENDDO
     492        ENDDO
    493493
    494494        DO k = 1, klev
    495495           DO i = 1, klon
    496               pi_var_day(i,k,id_aero) = pi_var_day(i,k,id_aero)*pi_load_src(i)/load_tgt(i)
    497            END DO
    498         END DO
     496              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))
     497           ENDDO
     498        ENDDO
    499499
    500500        IF (debug) THEN
    … …  
    504504                 zrho = pplay(i,k)/t_seri(i,k)/RD          ! [kg/m3]
    505505                 volm = pi_var_day(i,k,id_aero)*1.E-9/zrho ! [kg/kg]
    506                  load_tgt_test(i) = load_tgt_test(i) + 1/RG * volm * delp(i,k)
    507               END DO
    508            END DO
     506                 load_tgt_test(i) = load_tgt_test(i) + volm*delp(i,k)/RG
     507              ENDDO
     508           ENDDO
    509509           CALL writefield_phy('pi_day_tgt2',pi_var_day(:,:,id_aero),klev)
    510510           CALL writefield_phy('pi_load_tgt',load_tgt(:),1)
    511511           CALL writefield_phy('pi_load_tgt_test',load_tgt_test(:),1)
    512512           CALL writefield_phy('pi_load_src',pi_load_src(:),1)
    513         END IF
     513        ENDIF
    514514
    515515
    … …  
    519519        pi_var_day(:,:,id_aero) = tmp2(:,:)
    520520
    521      END IF ! vert_interp
     521     ENDIF ! vert_interp
    522522
    523523
    … …  
    539539                         trim(name_aero(id_aero)),'(i,k,im)=',           &
    540540                         var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero)
    541                  END IF
     541                 ENDIF
    542542                 WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero)
    543543                 WRITE(lunout,*) 'day1, day2, jDay = ', day1, day2, jDay
    544544                 CALL abort_physic('readaerosol_interp','Error in interpolation 1',1)
    545               END IF
    546            END DO
    547         END DO
    548      END IF
     545              ENDIF
     546           ENDDO
     547        ENDDO
     548     ENDIF
    549549
    550550     IF (MINVAL(pi_var_day(:,:,id_aero)) < 0. ) THEN
    … …  
    558558                         trim(name_aero(id_aero)),'(i,k,im)=',           &
    559559                         pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero)
    560                  END IF
     560                 ENDIF
    561561                 
    562562                 WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero)
    563563                 CALL abort_physic('readaerosol_interp','Error in interpolation 2',1)
    564               END IF
    565            END DO
    566         END DO
    567      END IF
    568 
    569   END IF ! lnewday
     564              ENDIF
     565           ENDDO
     566        ENDDO
     567     ENDIF
     568
     569  ENDIF ! lnewday
    570570
    571571!****************************************************************************************
Note: See TracChangeset for help on using the changeset viewer.