Changeset 2819

Timestamp:

Mar 8, 2017, 3:34:12 PM (8 years ago)

Author:

dcugnet

Message:

• Modification to allow previous convention climoz files (360 records for all calendars) to be handled again:
  • read_climoz=1 / 2 (tro3 only / tro3 and tro3_daylight) => 360 records file
  • read_climoz=3 / 4 (same) => daily (calendar compliant) file or monthly (14 records) file
• Few fixes in the ozone forcing reading routine

Location:

LMDZ5/trunk/libf/phylmd

Files:

• physiq_mod.F90 (modified) (3 diffs)
• regr_pr_time_av_m.F90 (modified) (12 diffs)

LMDZ5/trunk/libf/phylmd/physiq_mod.F90

@@ -1021 +1021 @@
                                   = ["tro3         ","tro3_daylight"]
     ! vars_climoz(1:read_climoz): variables names in climoz file.
+    ! vars_climoz(1:read_climoz-2) if read_climoz>2 (temporary)
+    REAL :: ro3i ! 0<=ro3i<=360 ; required time index in NetCDF file for
+                 ! the ozone fields, old method.
 
     include "YOMCST.h"
@@ -1972 +1975 @@
 
           wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
+       ELSE IF(ANY([1,2]==read_climoz)) THEN
+          ro3i = REAL(INT((days_elapsed+jh_cur-jh_1jan)/year_len*360.))
+          IF(INT(ro3i) == 360) ro3i = 359.
+          ! (This should never occur, except perhaps because of roundup
+          ! error. See documentation.)
+          WRITE(*,*)' >> Using 360 records ozone files (old method).'
+          CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),      &
+                               ro3i+1., press_edg_climoz, paprs, wo)
+          FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
+               * zmasse / dobson_u / 1e3
        ELSE
           IF(adjust_tropopause) THEN
@@ -1977 +1990 @@
             CALL dyn_tropopause(t_seri, ztsol, paprs, pplay, presnivs, rot,  &
                  pr_tropopause)
-            CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),    &
+            CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz-2),  &
                  days_elapsed+jh_cur-jh_1jan, press_edg_climoz, paprs, wo,   &
                  time_climoz,  latitude_deg,  press_cen_climoz, pr_tropopause)
           ELSE
             WRITE(*,*)' >> Interpolating O3 field directly on gcm levels.'
-            CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),    &
+            CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz-2),  &
                  days_elapsed+jh_cur-jh_1jan, press_edg_climoz, paprs, wo,   &
                  time_climoz)

LMDZ5/trunk/libf/phylmd/regr_pr_time_av_m.F90

@@ -71 +71 @@
   LOGICAL, SAVE :: lO3Trop                         !--- Tropopause ozone flag
   LOGICAL, SAVE :: lfirst=.TRUE.                   !--- First call flag
+!$OMP THREADPRIVATE(lfirst)
   REAL, PARAMETER :: co3(3)=[91.,28.,20.]          !--- Coeffs for o3 threshold
   REAL, PARAMETER :: prtrop=5.E+3  !--- Value lower than the tropopause pressure
@@ -79 +80 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE regr_pr_time_av(fID, nam, julien, pint_in, pcen_ou, v3, time_in,    &
+SUBROUTINE regr_pr_time_av(fID, nam, julien, pint_in, pint_ou, v3, time_in,    &
                                      lat_in, pcen_in, ptrop_ou)
 !
@@ -102 +103 @@
   CHARACTER(LEN=13), INTENT(IN) :: nam(:)     !--- NetCDF variables names
   REAL,    INTENT(IN)  :: julien        !--- Days since Jan 1st
-  REAL,    INTENT(IN)  :: pint_in(:)    !--- Interfaces Pressure, Pa, ascending
-  REAL,    INTENT(IN)  :: pcen_ou(:,:)  !--- Mid-layers LMDZ Pres, Pa (klon,llm+1)
+  REAL,    INTENT(IN)  :: pint_in(:)    !--- Interfaces file Pres, Pa, ascending
+  REAL,    INTENT(IN)  :: pint_ou(:,:)  !--- Interfaces LMDZ Pres, Pa (klon,llm+1)
   REAL,    INTENT(OUT) :: v3(:,:,:)     !--- Regridded field (klon,llm,SIZE(nam))
   REAL,    INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days
                                               !    since Jan 1 of current year
   REAL,    INTENT(IN), OPTIONAL :: lat_in(:)  !--- Input/output latitudes vector
-  REAL,    INTENT(IN), OPTIONAL :: pcen_in(:) !--- Mid-layers pressure
+  REAL,    INTENT(IN), OPTIONAL :: pcen_in(:) !--- Mid-layers file pressure
   REAL,    INTENT(IN), OPTIONAL :: ptrop_ou(:)!--- Output tropopause pressure (klon)
 !-------------------------------------------------------------------------------
@@ -136 +137 @@
 !-------------------------------------------------------------------------------
   sub="regr_pr_time_av"
+  nlev_in=SIZE(pint_in)-1
   CALL assert(SIZE(v3,1)==klon,                 TRIM(sub)//" v3 klon")
   CALL assert(SIZE(v3,2)==nbp_lev,              TRIM(sub)//" v3 nbp_lev")
   n_var = assert_eq(SIZE(nam), SIZE(v3,3),      TRIM(sub)//" v3 n_var")
-  CALL assert(SHAPE(pcen_ou)==[klon, nbp_lev+1],TRIM(sub)//" pcen_ou")
+  CALL assert(SHAPE(pint_ou)==[klon, nbp_lev+1],TRIM(sub)//" pint_ou")
   IF(PRESENT(lat_in))   CALL assert(SIZE(lat_in  )==klon,TRIM(sub)//" lat_in klon")
   IF(PRESENT(ptrop_ou)) CALL assert(SIZE(ptrop_ou)==klon,TRIM(sub)//" ptrop_ou klon")
+  IF(PRESENT(pcen_in))  CALL assert(SIZE(pcen_in)==nlev_in,TRIM(sub)//" pcen_in")
   ltrop=PRESENT(lat_in).AND.PRESENT(pcen_in).AND.PRESENT(ptrop_ou)
-  nlev_in=SIZE(pint_in)-1
 
   !$OMP MASTER
@@ -162 +164 @@
         IF(lO3Trop) ALLOCATE(otm(nbp_lon,nbp_lat),otp(nbp_lon,nbp_lat))
       END IF
-      IF(PRESENT(pcen_in)) itrp0=locate(pcen_in,prtrop) !--- Above tropopause: 50hPa
+      !--- INITIAL INDEX: LOCATE A LAYER WELL ABOVE TROPOPAUSE (50hPa)
+      IF(PRESENT(pcen_in)) THEN
+        itrp0=locate(pcen_in,prtrop)
+      ELSE
+        itrp0=locate(SQRT(pint_in(1:nlev_in)*pint_in(2:nlev_in+1)),prtrop)
+      END IF
       IF(lPrSurf) WRITE(*,*)' >> Using GROUND PRESSURE from input O3 forcing file.'
       IF(linterp) WRITE(*,*)' >> Monthly O3 files => ONLINE TIME INTERPOLATION.'
@@ -178 +185 @@
     !=== TIME INTERPOLATION FOR MONTHLY INPUT FILES
     IF(linterp) THEN
-      WRITE(*,'(3(a,f7.3))')'  >> Interpolating O3 at julian day ',julien,' fr&
-      &om forcing fields at times ',time_in(irec),' and ', time_in(irec+1)
+      WRITE(*,'(3(a,f7.3))')'  >> Interpolating O3 at julian day ',julien,     &
+      &' from forcing fields at times ',time_in(irec),' and ', time_in(irec+1)
       al=(time_in(irec+1)-julien)/(time_in(irec+1)-time_in(irec))
       v1=al*v1m+(1.-al)*v1p
@@ -188 +195 @@
   END IF
   !$OMP END MASTER
-  !$OMP BARRIER
   CALL scatter2d(v1,v2)
-  IF(PRESENT(pcen_in)) CALL bcast(itrp0)
-  !--- "ps" not in input file => assume it is equal to current LMDZ value
-  IF(lPrSurf) THEN; CALL scatter2d(ps1,ps2); ELSE; ps2=pcen_ou(:,1); END IF
+  CALL bcast(itrp0)
+  !--- No "ps" in input file => assumed to be equal to current LMDZ ground press
+  IF(lPrSurf) THEN; CALL scatter2d(ps1,ps2); ELSE; ps2=Pint_ou(:,1); END IF
   IF(lPrTrop) CALL scatter2d(pt1,pt2)
   IF(lO3Trop) CALL scatter2d(ot1,ot2)
@@ -199 +205 @@
   IF(.NOT.ltrop) THEN
     DO i=1,klon
-      CALL regr_conserv(1,v2(i,:,:) , Pint_in , Pcen_ou(i,nbp_lev+1:1:-1),    &
+      CALL regr_conserv(1,v2(i,:,:) , Pint_in , Pint_ou(i,nbp_lev+1:1:-1),    &
                           v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),pint_in))
     END DO
@@ -206 +212 @@
     DO i=1,klon
       SigT_in = get_SigTrop(i,itrp)        !--- input  (file) tropopause
-      SigT_ou = ptrop_ou(i)/pcen_ou(i,1)   !--- output (lmdz) tropopause
+      SigT_ou = ptrop_ou(i)/pint_ou(i,1)   !--- output (lmdz) tropopause
       SigT_mn = SQRT(SigT_in*SigT_ou)      !--- mean tropopause>strained domain
 
@@ -221 +227 @@
       !--- STRAINED INPUT logP PROFILE ; alpha: MAX. STRETCH. EXPONENT INCREMENT
       alpha = LOG(SigT_ou/SigT_in)/LOG(SigT_in)
-      Pint_st(:) = pcen_ou(i,1) * Sig_in(:)**(1+alpha*phi(:))
+      Pint_st(:) = pint_ou(i,1) * Sig_in(:)**(1+alpha*phi(:))
 
       !--- REGRID INPUT PROFILE ON STRAINED VERTICAL LEVELS
-      CALL regr_conserv(1,v2(i,:,:) , Pint_st,  Pcen_ou(i,nbp_lev+1:1:-1),     &
+      CALL regr_conserv(1,v2(i,:,:) , Pint_st,  Pint_ou(i,nbp_lev+1:1:-1),     &
                           v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),Pint_st))
     END DO
@@ -275 +281 @@
   END IF
   lfirst=.FALSE.
+  CALL bcast(lfirst)
 
 END SUBROUTINE update_fields
@@ -361 +368 @@
   REAL    :: prt                               !--- Air pressure   at tropopause
   REAL    :: al                                !--- Interpolation coefficient
-  REAL    :: y_frac                            !--- Elapsed year fraction
   REAL    :: coef                              !--- Coef: North/South transition
 !-------------------------------------------------------------------------------