Changeset 1425 for LMDZ4/trunk/libf/dyn3d/startvar.F90

Timestamp:

Sep 2, 2010, 3:45:23 PM (14 years ago)

Author:

lguez

Message:

Replaced Numerical Recipes procedures for spline interpolation (not in
the public domain) by procedures from the Pchip package in the Slatec
library. This only affects the program "ce0l", not the program
"gcm". Tested on Brodie SX8 with "-debug" and "-prod", "-parallel
none" and "-parallel mpi". "start.nc" and "limit.nc" are
changed. "startphy.nc" is not changed. The relative change is of order
1e-7 or less. The revision makes the program faster (tested on Brodie
with "-prod -d 144x142x39", CPU time is 38 s, instead of 54
s). Procedures from Slatec are untouched, except for
"i1mach.F". Created procedures "pchfe_95" and "pchsp_95" which are
wrappers for "pchfe" and "pchsp" from Slatec. "pchfe_95" and
"pchsp_95" have a safer and simpler interface.

Replaced "make" by "sxgmake" in "arch-SX8_BRODIE.fcm". Added files for
compilation by FCM with "g95".

In "arch-linux-32bit.fcm", replaced "pgf90" by "pgf95". There was no
difference between "dev" and "debug" so added "-O1" to "dev". Added
debugging options. Removed "-Wl,-Bstatic
-L/usr/lib/gcc-lib/i386-linux/2.95.2", which usually produces an error
at link-time.

Bash is now ubiquitous while KornShell? is not so use Bash instead of
KornShell? in FCM.

Replaced some statements "write(6,*)" by "write(lunout,*)". Replaced
"stop" by "stop 1" in the case where "abort_gcm" is called with "ierr
/= 0". Removed "stop" statements at the end of procedures
"limit_netcdf" and main program "ce0l" (why not let the program end
normally?).

Made some arrays automatic instead of allocatable in "start_inter_3d".

Zeroed "wake_pe", "fm_therm", "entr_therm" and "detr_therm" in
"dyn3dpar/etat0_netcdf.F90". The parallel and sequential results of
"ce0l" are thus identical.

File:

• LMDZ4/trunk/libf/dyn3d/startvar.F90 (modified) (2 diffs)

LMDZ4/trunk/libf/dyn3d/startvar.F90

@@ -643 +643 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2,lon_in2,&
-                          lat_in2, pls_in, var3d, ibar)
-!
-!-------------------------------------------------------------------------------
+SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2, &
+     lon_in2, lat_in2, pls_in, var3d, ibar)
+
+  use pchsp_95_m, only: pchsp_95
+  use pchfe_95_m, only: pchfe_95
+
 ! Arguments:
   CHARACTER(LEN=*),             INTENT(IN)    :: varname
@@ -655 +657 @@
   REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in2
   REAL, DIMENSION(jml2),        INTENT(IN)    :: lat_in2
-  REAL, DIMENSION(iml,jml,lml), INTENT(IN)    :: pls_in
-  REAL, DIMENSION(iml,jml,lml), INTENT(OUT)   :: var3d
+  REAL, DIMENSION(iml, jml, lml), INTENT(IN)    :: pls_in
+  REAL, DIMENSION(iml, jml, lml), INTENT(OUT)   :: var3d
   LOGICAL,                      INTENT(IN)    :: ibar
-!-------------------------------------------------------------------------------
+!----------------------------------------------------------------------------
 ! Local variables:
 #include "iniprint.h"
-  LOGICAL               :: check=.TRUE.
-  REAL                  :: bx, by, chmin, chmax
-  INTEGER               :: ii, ij, il
-  REAL,    DIMENSION(:,:,:), ALLOCATABLE :: var_tmp3d
+  LOGICAL:: check=.TRUE., skip
+  REAL                  chmin, chmax
+  INTEGER ii, ij, il, ierr
+  integer n_extrap ! number of extrapolated points
+  REAL, DIMENSION(iml, jml, llm_dyn):: var_tmp3d
   REAL,    DIMENSION(:),     ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
-  REAL,    DIMENSION(:),     ALLOCATABLE :: lev_dyn, ax, ay, yder
-  INTEGER, DIMENSION(:),     ALLOCATABLE :: lind
-!-------------------------------------------------------------------------------
-  IF(check) WRITE(lunout,*)'Going into flinget to extract the 3D  field.'
-  IF(check) WRITE(lunout,*)fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn
-  IF(check) WRITE(lunout,*)'Allocating space for interpolation',iml,jml,llm_dyn
-
-  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn,jml_dyn,llm_dyn))
-  CALL flinget(fid_dyn,varname,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d)
+  REAL, DIMENSION(llm_dyn):: lev_dyn, ax, ay, yder
+
+!---------------------------------------------------------------------------
+  IF(check) WRITE(lunout, *)'Going into flinget to extract the 3D  field.'
+  IF(check) WRITE(lunout, *) fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, &
+       ttm_dyn
+  IF(check) WRITE(lunout, *) 'Allocating space for interpolation', iml, jml, &
+       llm_dyn
+
+  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
+  CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, 1, 1, &
+       var_ana3d)
 
 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
-  ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
-  lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
-  lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
+  ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn))
+  lon_ini(:)=lon_dyn(:, 1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
+  lat_ini(:)=lat_dyn(1, :); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
 
 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
-  ALLOCATE(lon_rad(iml_dyn),lat_rad(jml_dyn),lev_dyn(llm_dyn))
-  CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini,       &
+  ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn))
+  CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini,      &
                    levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
-  DEALLOCATE(lon_ini,lat_ini)
+  DEALLOCATE(lon_ini, lat_ini)
 
 !--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
-  ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
-  DO il=1,llm_dyn
-    CALL interp_startvar(varname, ibar, il==1,                                 &
-      iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana3d(:,:,il), iml, jml, jml2,   &
-      lon_in,  lat_in,  lon_in2, lat_in2, var_tmp3d(:,:,il))
+  DO il=1, llm_dyn
+    CALL interp_startvar(varname, ibar, il==1, iml_dyn, jml_dyn, lon_rad, &
+         lat_rad, var_ana3d(:, :, il), iml, jml, jml2, lon_in, lat_in, &
+         lon_in2, lat_in2, var_tmp3d(:, :, il))
   END DO
-  DEALLOCATE(lon_rad,lat_rad)
-
-  ALLOCATE(lind(llm_dyn))
-  DO il=1,llm_dyn
-    lind(il) = llm_dyn-il+1
+  DEALLOCATE(lon_rad, lat_rad)
+
+!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
+  ax = lev_dyn(llm_dyn:1:-1) 
+  skip = .false.
+  n_extrap = 0
+  DO ij=1, jml
+    DO ii=1, iml-1
+      ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
+      yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.)
+      CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), &
+           var3d(ii, ij, lml:1:-1), ierr)
+      if (ierr < 0) stop 1
+      n_extrap = n_extrap + ierr
+    END DO
   END DO
-
-!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
-  ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn))
-  DO ij=1,jml
-    DO ii=1,iml-1
-      ax(:)=lev_dyn(lind(:)) 
-      ay(:)=var_tmp3d(ii,ij,lind(:))
-      CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder)
-      DO il=1,lml
-        bx=pls_in(ii,ij,il)
-        CALL SPLINT(ax, ay, yder, llm_dyn, bx, by)
-        var3d(ii,ij,il) = by
-      END DO
-    END DO
-    var3d(iml,ij,:) = var3d(1,ij,:) 
+  if (n_extrap /= 0) then
+     print *, "start_inter_3d pchfe_95: n_extrap = ", n_extrap
+  end if
+  var3d(iml, :, :) = var3d(1, :, :) 
+
+  DO il=1, lml
+    CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
+    WRITE(lunout, *)' '//TRIM(varname)//'  min max l ', il, chmin, chmax
   END DO
-  DEALLOCATE(var_tmp3d,lev_dyn,ax,ay,yder,lind)
-
-  DO il=1,lml
-    CALL minmax(iml*jml,var3d(1,1,il),chmin,chmax)
-    WRITE(lunout,*)' '//TRIM(varname)//'  min max l ',il,chmin,chmax
-  END DO
-
-  RETURN
 
 END SUBROUTINE start_inter_3d