Changeset 2542 for LMDZ5/branches/testing/libf/dynphy_lonlat

Timestamp:

Jun 6, 2016, 4:04:57 PM (8 years ago)

Author:

Laurent Fairhead

Message:

Merged trunk changes r2487:2541 into testing branch

Location:

LMDZ5/branches/testing

Files:

• . (modified) (1 prop)
• libf/dynphy_lonlat/grid_atob_m.f90 (modified) (3 diffs)
• libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90 (modified) (1 diff)
• libf/dynphy_lonlat/phylmd/limit_netcdf.F90 (modified) (9 diffs)

LMDZ5/branches/testing/libf/dynphy_lonlat/grid_atob_m.f90

@@ -29 +29 @@
 ! Local variables:
   CHARACTER(LEN=256) :: modname="fine2coarse"
-  INTEGER :: mi, ni, ii, ji, mo, no, io, jo, nr(2), m1, n1, m2, n2, mx, my, nn
-  LOGICAL :: li, lo
-  REAL    :: inc
+  INTEGER :: mi, ni, ii, ji, mo, no, io, jo, nr(2), m1,m2, n1,n2, mx,my, nn, i,j
+  LOGICAL :: li, lo, first=.TRUE.
+  REAL    :: inc, cpa, spa, crlo(SIZE(x_i))
+  REAL, SAVE :: pi, hpi
   INTEGER, DIMENSION(SIZE(x_o),SIZE(y_o)) :: num_tot
   LOGICAL, DIMENSION(SIZE(x_o),SIZE(y_o)) :: found, mask
-  REAL,    DIMENSION(SIZE(x_o),SIZE(y_o)) :: dist
+  REAL,    DIMENSION(SIZE(x_i),SIZE(y_i)) :: dist
   REAL,    DIMENSION(SIZE(x_o))           :: a, b
   REAL,    DIMENSION(SIZE(y_o))           :: c, d
   REAL,    PARAMETER :: thresh=1.E-5
 !-------------------------------------------------------------------------------
+  IF(first) THEN; pi=4.0*ATAN(1.0); hpi=pi/2.0; first=.FALSE.; END IF
   mi=SIZE(x_i); ni=SIZE(y_i); mo=SIZE(x_o); no=SIZE(y_o)
   m1=m1; m2=mo; mx=mo; IF(PRESENT(msk)) mx=SIZE(msk,1)
@@ -91 +93 @@
       IF(found(io,jo)) CYCLE
 !      IF(prt_level>=1) PRINT*, "Problem: point out of domain (i,j)=", io,jo
-      CALL dist_sphe(x_o(io),y_o(jo),x_i,y_i,dist(:,:))
+      crlo(:)=COS(x_o(io)-x_i(:))     !--- COS of points 1 and 2 angle
+      cpa=COS(y_o(jo)); spa=SIN(y_o(jo))
+      DO j=1,ni; dist(:,j)=ACOS(spa*SIN(y_i(j))+cpa*COS(y_i(j))*crlo(:)); END DO
       nr=MINLOC(dist(:,:))!; IF(prt_level>=1) PRINT*, "Solution: ", nr
       inc=1.0; IF(li) inc=d_i(nr(1),nr(2))
@@ -252 +256 @@
 !-------------------------------------------------------------------------------
 
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE dist_sphe(rf_lon,rf_lat,rlon,rlat,distance)
-!
-!-------------------------------------------------------------------------------
-! Author:  Laurent Li (december 30th 1996).
-! Purpose: Compute min. distance (along big circle) between 2 points in radians.
-!-------------------------------------------------------------------------------
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  REAL, INTENT(IN) :: rf_lon, rf_lat  !--- Reference point coordinates (radians)
-  REAL, INTENT(IN) :: rlon(:), rlat(:)!--- Points longitudes/latitudes (radians)
-  REAL, INTENT(OUT):: distance(SIZE(rlon),SIZE(rlat)) !--- Distance    (radians)
-!-------------------------------------------------------------------------------
-! Local variables:
-  LOGICAL, SAVE :: first=.TRUE.
-  REAL,    SAVE :: pi, hpi
-  REAL    :: pa, pb, cpa, spa, crlo(SIZE(rlon))
-  INTEGER :: i, j
-!-------------------------------------------------------------------------------
-  IF(first) THEN; pi=4.0*ATAN(1.0); hpi=pi/2.0; first=.FALSE.; END IF
-  crlo(:)=COS(rf_lon-rlon(:))         !--- COS of points 1 and 2 angle
-  pa=hpi-rf_lat                       !--- North Pole - Point 1 distance
-  cpa=COS(pa); spa=SIN(pa)
-  DO j=1,SIZE(rlat)
-    pb=hpi-rlat(j)                    !--- North Pole - Point 2 distance
-    distance(:,j)=ACOS(cpa*COS(pb)+spa*SIN(pb)*crlo(:))
-  END DO
-
-END SUBROUTINE dist_sphe
-!
-!-------------------------------------------------------------------------------
-
 END MODULE grid_atob_m
 !

LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90

@@ -111 +111 @@
   LOGICAL :: ok_ade, ok_aie, ok_cdnc, aerosol_couple
   INTEGER :: flag_aerosol
-  LOGICAL :: flag_aerosol_strat
+  INTEGER :: flag_aerosol_strat
   LOGICAL :: new_aod
   REAL    :: bl95_b0, bl95_b1

LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/limit_netcdf.F90

@@ -60 +60 @@
 !  *    11/2009: L. Guez     (ozone day & night climatos, see etat0_netcdf.F90)
 !  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
+!  *    04/2016: D. Cugnet   (12/14 recs SST/SIC files: cyclic/interannual runs)
 !-------------------------------------------------------------------------------
 #ifndef CPP_1D
@@ -339 +340 @@
   INTEGER           :: ndays_in           ! number of days
 !--- misc
-  INTEGER           :: i, j, k, l         ! loop counters
+  INTEGER           :: i, j, k, l, ll     ! loop counters
   REAL, ALLOCATABLE :: work(:,:)          ! used for extrapolation
-  CHARACTER(LEN=25) :: title              ! for messages
+  CHARACTER(LEN=128):: title, mess        ! for messages
   LOGICAL           :: extrp              ! flag for extrapolation
   REAL              :: chmin, chmax
@@ -392 +393 @@
 !--- Time (variable is not needed - it is rebuilt - but calendar is)
   CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep), fnam)
-  ALLOCATE(timeyear(lmdep))
+  ALLOCATE(timeyear(MAX(lmdep,14)))
   CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
   cal_in=' '
@@ -405 +406 @@
   CALL msg(5,'var, calendar, dim: '//TRIM(dnam)//' '//TRIM(cal_in), lmdep)
   
-!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION --------------------
-  !--- Determining input file number of days, depending on calendar
+!--- Determining input file number of days, depending on calendar
   ndays_in=year_len(anneeref, cal_in)
 
-!--- Time vector reconstruction (time vector from file is not trusted)
-!--- If input records are not monthly, time sampling has to be constant !
-  timeyear=mid_months(anneeref, cal_in, lmdep)
-  IF (lmdep /= 12) WRITE(lunout,*) 'Note : les fichiers de ', TRIM(mode), &
-       ' ne comportent pas 12, mais ', lmdep, ' enregistrements.'
+!--- Rebuilding input time vector (field from input file might be unreliable)
+  IF(lmdep==12) THEN
+    timeyear=mid_month(anneeref, cal_in, ndays_in)
+    CALL msg(0,'Monthly periodic input file (perpetual run).')
+  ELSE IF(lmdep==14) THEN
+    timeyear=mid_month(anneeref, cal_in, ndays_in)
+    CALL msg(0,'Monthly 14-records input file (interannual run).')
+  ELSE IF(lmdep==ndays_in) THEN
+    timeyear=[(REAL(k)-0.5,k=1,ndays_in)]
+    CALL msg(0,'Daily input file (no time interpolation).')
+  ELSE
+    WRITE(mess,'(a,i3,a,i3,a)')'Mismatching input file: found',lmdep,        &
+      ' records, 12/14/',ndays_in,' (periodic/interannual/daily) needed'
+    CALL abort_physic('mid_months',TRIM(mess),1)
+  END IF
 
 !--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
-  ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep))
+  ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, MAX(lmdep,14)))
   IF(extrp) ALLOCATE(work(imdep, jmdep))
   CALL msg(5,'')
@@ -436 +446 @@
       WHERE(NINT(mask)/=1) champint=0.001
     END IF
-    champtime(:, :, l)=champint
+    !--- Special case for periodic input file: index shifted
+    ll = l; IF(lmdep==12) ll = l + 1
+    champtime(:, :, ll)=champint
   END DO
+  IF(lmdep==12) THEN
+    champtime(:,:, 1)=champtime(:,:,13)
+    champtime(:,:,14)=champtime(:,:, 2)
+  END IF
   CALL ncerr(NF90_CLOSE(ncid), fnam)
 
@@ -444 +460 @@
 
 !--- TIME INTERPOLATION ------------------------------------------------------
-  IF(prt_level>5) THEN
-     WRITE(lunout, *)
-     WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
-     WRITE(lunout, *)' Vecteur temps en entree: ', timeyear
-     WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays
-  END IF
-
-  ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays))
-  skip = .false.
-  n_extrap = 0
-  DO j=1, jjp1
-    DO i=1, iim
-      yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, &
-           vc_beg=0., vc_end=0.)
-      CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
-           arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
-      if (ierr < 0) stop 1
-      n_extrap = n_extrap + ierr
-    END DO
-  END DO
-  if (n_extrap /= 0) then
-     WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
-  end if
+  IF(prt_level>0) THEN
+     IF(ndays/=ndays_in) THEN
+        WRITE(lunout, *)
+        WRITE(lunout,*)'DIFFERENTES LONGEURS D ANNEES:'
+        WRITE(lunout,*)' Dans le fichier d entree: ',ndays_in
+        WRITE(lunout,*)' Dans le fichier de sortie: ',ndays
+     END IF
+     IF(lmdep==ndays_in) THEN
+        WRITE(lunout, *)
+        WRITE(lunout, *)'PAS D INTERPOLATION TEMPORELLE.'
+        WRITE(lunout, *)' Fichier journalier en entree.'
+     ELSE
+        WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
+        WRITE(lunout, *)' Vecteur temps en entree: ', timeyear
+        WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays-1
+     END IF
+  END IF
+  ALLOCATE(yder(MAX(lmdep,14)), champan(iip1, jjp1, ndays))
+  IF(lmdep==ndays_in) THEN
+     champan(1:iim,:,:)=champtime
+  ELSE
+     skip = .false.
+     n_extrap = 0
+     DO j=1, jjp1
+       DO i=1, iim
+         yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, &
+              vc_beg=0., vc_end=0.)
+         CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
+              arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
+         if (ierr < 0) stop 1
+         n_extrap = n_extrap + ierr
+       END DO
+     END DO
+     IF(n_extrap /= 0) WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
+  END IF
   champan(iip1, :, :)=champan(1, :, :)
   DEALLOCATE(yder, champtime, timeyear)
@@ -752 +780 @@
 !-------------------------------------------------------------------------------
 !
-FUNCTION mid_months(y,cal_in,nm)
+FUNCTION mid_month(y,cal_in,ndays_out)
 !
 !-------------------------------------------------------------------------------
@@ -758 +786 @@
 !-------------------------------------------------------------------------------
 ! Arguments:
-  INTEGER,                INTENT(IN) :: y               ! year
-  CHARACTER(LEN=*),       INTENT(IN) :: cal_in          ! calendar
-  INTEGER,                INTENT(IN) :: nm              ! months/year number
-  REAL,    DIMENSION(nm)             :: mid_months      ! mid-month times
+  INTEGER,          INTENT(IN) :: y           ! year
+  CHARACTER(LEN=*), INTENT(IN) :: cal_in      ! calendar
+  INTEGER,          INTENT(IN) :: ndays_out   ! days   number
+  REAL,          DIMENSION(14) :: mid_month   ! mid-bins times
 !-------------------------------------------------------------------------------
 ! Local variables:
-  CHARACTER(LEN=99)                  :: mess            ! error message
-  CHARACTER(LEN=20)                  :: cal_out         ! calendar (for outputs)
-  INTEGER, DIMENSION(nm)             :: mnth            ! months lengths (days)
-  INTEGER                            :: m               ! months counter
-  INTEGER                            :: nd              ! number of days
-!-------------------------------------------------------------------------------
-  nd=year_len(y,cal_in)
-
-  IF(nm==12) THEN
-
-  !--- Getting the input calendar to reset at the end of the function
-    CALL ioget_calendar(cal_out)
-
-  !--- Unlocking calendar and setting it to wanted one
-    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
-
-  !--- Getting the length of each month
-    DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO
+  CHARACTER(LEN=99)      :: mess              ! error message
+  CHARACTER(LEN=20)      :: cal_out           ! output calendar
+  INTEGER, DIMENSION(14) :: tlen              ! months lengths (days)
+  INTEGER                :: m                 ! months counter
+  INTEGER                :: nd                ! number of days
+!-------------------------------------------------------------------------------
+  !--- Save the input calendar to reset it at the end of the function
+  CALL ioget_calendar(cal_out)
+
+  !--- Unlock calendar and set it to "cal_in"
+  CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
+
+  !--- Get the length of each month
+  tlen(1 )=ioget_mon_len(y-1,12)
+  DO m=1,12; tlen(m+1)=ioget_mon_len(y,m); END DO
+  tlen(14)=ioget_mon_len(y+1, 1)
+
+  !--- Mid-bins times
+  mid_month(1)=-0.5*REAL(tlen(1))
+  DO m=2,14; mid_month(m)=mid_month(m-1)+0.5*REAL(tlen(m-1)+tlen(m)); END DO
 
   !--- Back to original calendar
-    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
-
-  ELSE IF(MODULO(nd,nm)/=0) THEN
-    WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',&
-      nm,' months/year. Months number should divide days number.'
-    CALL abort_physic('mid_months',TRIM(mess),1)
-
-  ELSE
-    mnth=[(m,m=1,nm,nd/nm)]
-  END IF
-
-!--- Mid-months times
-  mid_months(1)=0.5*REAL(mnth(1))
-  DO k=2,nm
-    mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k))
-  END DO
-
-END FUNCTION mid_months
+  CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
+
+END FUNCTION mid_month
 !
 !-------------------------------------------------------------------------------
@@ -863 +878 @@
 !
 !*******************************************************************************
+