Changeset 2820 for LMDZ5/trunk/libf/phylmd/regr_pr_time_av_m.F90

Timestamp:

Mar 10, 2017, 5:46:08 PM (8 years ago)

Author:

dcugnet

Message:

Fix and cleaning/improvements over the previous commit:

• Old style climoz files are handled without trouble.
• New style climoz (monthly with 14 records or daily) are handled

with the same read_climoz key values (1 or 2).

• Daily files reading activation is done using the same key

used in ce0l to create them: ok_daily_climoz=y. This key is now
put to FALSE by default (old style 360 days files are not considered
as daily).

File:

• LMDZ5/trunk/libf/phylmd/regr_pr_time_av_m.F90 (modified) (13 diffs)

LMDZ5/trunk/libf/phylmd/regr_pr_time_av_m.F90

@@ -2 +2 @@
 MODULE regr_pr_time_av_m
 
+  USE print_control_mod, ONLY: lunout
   USE write_field_phy
   USE mod_phys_lmdz_transfert_para, ONLY: bcast
@@ -117 +118 @@
   CHARACTER(LEN=80)  :: sub
   CHARACTER(LEN=320) :: msg
-  INTEGER :: vID, ncerr, n_var, nlev_in, ndim, i, ibot, itop, itrp, itrp0
-  LOGICAL :: ltrop                            !--- Need to adjust tropopause
+  INTEGER :: vID, ncerr, n_var, nlev_in,ntim_in, ndim, i, ibot, itop, itrp,itrp0
+  LOGICAL :: lAdjTro                          !--- Need to adjust tropopause
   REAL    :: y_frac                           !--- Elapsed year fraction
   REAL    :: alpha, beta, al                  !--- For strectching/interpolation
@@ -137 +138 @@
 !-------------------------------------------------------------------------------
   sub="regr_pr_time_av"
-  nlev_in=SIZE(pint_in)-1
+  nlev_in=SIZE(pint_in)-1; ntim_in=SIZE(time_in)
   CALL assert(SIZE(v3,1)==klon,                 TRIM(sub)//" v3 klon")
   CALL assert(SIZE(v3,2)==nbp_lev,              TRIM(sub)//" v3 nbp_lev")
@@ -145 +146 @@
   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)
+  lAdjTro=PRESENT(ptrop_ou)
+  IF(lAdjTro.AND.(.NOT.PRESENT(lat_in).OR..NOT.PRESENT(pcen_in))) &
+  CALL abort_physic(sub, 'Missing lat_in and/or pcen_in (adjust_tropopause=T)', 1)
 
   !$OMP MASTER
@@ -155 +158 @@
       lPrTrop=NF90_INQ_VARID(fID,"tropopause_air_pressure",vID)==NF90_NOERR
       lO3Trop=NF90_INQ_VARID(fID,"tro3_at_tropopause"     ,vID)==NF90_NOERR
-      linterp=PRESENT(time_in)
-      IF(linterp) linterp=SIZE(time_in,1)==14
+      linterp=PRESENT(time_in); IF(linterp) linterp=ntim_in==14
       IF(linterp) THEN
         ALLOCATE(v1m(nbp_lon,nbp_lat,nlev_in,n_var))
@@ -165 +167 @@
       END IF
       !--- 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.'
-      WRITE(*,*)' >> o3 forcing file tropopause location uses:'
-      IF(lPrTrop) THEN;      WRITE(*,*)'    PRESSURE AT TROPOPAUSE from file.'
-      ELSE IF(lO3Trop) THEN; WRITE(*,*)'    O3 CONCENTRATION AT TROPOPAUSE from file.'
-      ELSE;                  WRITE(*,*)'    PARAMETRIZATION of O3 concentration at tropopause.'
+      IF(lAdjTro) itrp0=locate(pcen_in,prtrop)
+      IF(lPrSurf) WRITE(lunout,*)TRIM(sub)//': Using GROUND PRESSURE from input O3 forcing file.'
+      IF(linterp) WRITE(lunout,*)TRIM(sub)//': Monthly O3 files => ONLINE TIME INTERPOLATION.'
+      IF(lAdjTro) THEN
+        WRITE(lunout,*)TRIM(sub)//': o3 forcing file tropopause location uses:'
+        IF(lPrTrop) THEN
+          WRITE(lunout,*)'    PRESSURE AT TROPOPAUSE from file.'
+        ELSE IF(lO3Trop) THEN
+          WRITE(lunout,*)'    O3 CONCENTRATION AT TROPOPAUSE from file.'
+        ELSE
+          WRITE(lunout,*)'    PARAMETRIZED O3 concentration at tropopause.'
+        END IF
       END IF
     END IF
 
     !=== UPDATE (ALWAYS FOR DAILY FILES, EACH MONTH FOR MONTHLY FILES)
-    IF(.NOT.linterp.OR.(linterp.AND.(lfirst.OR.julien>=time_in(irec+1)))) &
-      CALL update_fields()
+    CALL update_fields()
+
 
     !=== TIME INTERPOLATION FOR MONTHLY INPUT FILES
     IF(linterp) THEN
-      WRITE(*,'(3(a,f7.3))')'  >> Interpolating O3 at julian day ',julien,     &
-      &' from forcing fields at times ',time_in(irec),' and ', time_in(irec+1)
+      WRITE(lunout,'(3(a,f7.3))')TRIM(sub)//': Interpolating O3 at julian day '&
+      &,julien,' from 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
@@ -195 +198 @@
   END IF
   !$OMP END MASTER
+  IF(lfirst) THEN
+    lfirst=.FALSE. ;     CALL bcast(lfirst) ; CALL bcast(lPrSurf)
+    CALL bcast(lPrTrop); CALL bcast(lO3Trop); CALL bcast(linterp)
+    IF(lAdjTro) CALL bcast(itrp0)
+  END IF
   CALL scatter2d(v1,v2)
-  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
@@ -203 +210 @@
 
   !--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing
-  IF(.NOT.ltrop) THEN
+  IF(.NOT.lAdjTro) THEN
     DO i=1,klon
       CALL regr_conserv(1,v2(i,:,:) , Pint_in , Pint_ou(i,nbp_lev+1:1:-1),    &
@@ -245 +252 @@
 !-------------------------------------------------------------------------------
   IF(.NOT.linterp) THEN                 !=== DAILY FILES: NO TIME INTERPOLATION
-    WRITE(*,*)' >> Updating Ozone forcing field: read from file.'
-    irec=INT(julien)+1                  !--- irec is just the julian day number
+    WRITE(lunout,*)TRIM(sub)//': Updating Ozone forcing field: read from file.'
+    irec=MIN(INT(julien)+1,ntim_in)     !--- irec is just the julian day number
+    !--- MIN -> Security in the unlikely case of roundup errors.
     CALL get_3Dfields(v1)               !--- Read ozone field(s)
-    IF(ltrop) THEN                      !--- Additional files for fields strain
+    IF(lAdjTro) THEN                    !--- Additional files for fields strain
       IF(lPrSurf) CALL get_2Dfield(ps1,"ps")
       IF(lPrTrop) CALL get_2Dfield(pt1,"tropopause_air_pressure")
@@ -254 +262 @@
     END IF
   ELSE                                  !=== MONTHLY FILES: GET 2 NEAREST RECS
-    WRITE(*,*)' >> Refreshing adjacent Ozone forcing fields.'
+    IF(.NOT.lfirst.AND.julien<time_in(irec+1)) RETURN
+    WRITE(lunout,*)TRIM(sub)//': Refreshing adjacent Ozone forcing fields.'
     IF(lfirst) THEN                     !=== READ EARLIEST TIME FIELDS
       irec=locate(time_in,julien)       !--- Need to locate surrounding times
       CALL get_3Dfields(v1m)            !--- Read ozone field(s)
-      IF(ltrop) THEN                    !--- Additional files for fields strain
+      IF(lAdjTro) THEN                  !--- Additional files for fields strain
         IF(lPrSurf) CALL get_2Dfield(psm,"ps")
         IF(lPrTrop) CALL get_2Dfield(ptm,"tropopause_air_pressure")
@@ -265 +274 @@
     ELSE                                !=== SHIFT FIELDS
       v1m=v1p                           !--- Ozone fields
-      IF(ltrop) THEN                    !--- Additional files for fields strain
+      IF(lAdjTro) THEN                  !--- Additional files for fields strain
         IF(lPrSurf) psm=psp             !--- Surface pressure
         IF(lPrTrop) ptm=ptp             !--- Tropopause pressure
@@ -273 +282 @@
     irec=irec+1
     CALL get_3Dfields(v1p)              !--- Read ozone field(s)
-    IF(ltrop) THEN                      !--- Additional files for fields strain
+    IF(lAdjTro) THEN                    !--- Additional files for fields strain
       IF(lPrSurf) CALL get_2Dfield(psp,"ps")
       IF(lPrTrop) CALL get_2Dfield(ptp,"tropopause_air_pressure")
@@ -280 +289 @@
     IF(lfirst) irec=irec-1
   END IF
-  lfirst=.FALSE.
-  CALL bcast(lfirst)
 
 END SUBROUTINE update_fields