!copyright (C) 2001 MSC-RPN COMM %%%RPNPHY%%%
***S/P INISURF
**
*
#include "phy_macros_f.h"
SUBROUTINE INISURF( E, ESIZ, F, FSIZ, NI, NK) 2,6
*
#include "impnone.cdk"
*
INTEGER NI, NK
INTEGER ESIZ, FSIZ
REAL E(ESIZ), F(FSIZ)
*
*
*Author
* Stephane Belair (February 1999)
*
*
*Revision
* 001 S. Belair (Mar 1999)
* New variable names
*
* 002 S. Belair (May 1999)
* Special treatment for the snow in ISBA
* 003 N. Brunet (Jul 1999)
* add several treatments on geophysical fields
* 004 S. Belair (Sept. 1999)
* New subroutine for coherence (soil and vegetation
* fields)
* 005 S. Belair (Feb 2000)
* Fix bug for the coherence between the soil texture
* and the mask MG.
* Sea ice and glaciers temperature is also affected by
* the filtering of mountains
* Code clean up
* 006 B. Bilodeau (Nov 2000)
* Change units of factor prcor (multiply by "grav")
* 007 B. Bilodeau (Jan 2001)
* New comdeck phybus.cdk. Change dimensions of a few
* arrays in the entry bus. Automatic arrays.
* 008 B. Bilodeau and S. Belair (May 2000) - Adaptation for vegetation
* in the southern hemisphere
* 009 B. Dugas (June 2000) - Initialize ALNOSNO for Force-Restore
* and modify call to CALCALB. Modify TGLACIER coherence test.
* 010 D. Talbot and B. Bilodeau (Oct 2001) - Add DHDX, DHDY and DHDXDY
*
*
*Object
* Transfer and initialize geophysical fields for the
* surface schemes
*
*
*Arguments
*
* - Input/Ouput -
* F field for permanent physics variables
* FSIZ dimension of F
* E field for entry variables
* ESIZ dimension of E
* NI horizontal dimension
*
*
**************************************************************
*NOTE: ***** This subroutine expects snow depth in cm.
* The snow depth is converted in metre (in this s/r)
* when the 'entry variables' are transfered to the
* permanent variables.
***************************************************************
*
#include "nclassvg.cdk"
*
**
#include "indx_sfc.cdk"
#include "consphy.cdk"
#include "isbapar.cdk"
#include "surfacepar.cdk"
*
#include "options.cdk"
*
#include "phybus.cdk"
#include "himin.cdk"
SAVE HIMIN
*
*
* the initial value of rugosity
* over water and ice
*
REAL Z0ICE, Z0SEA
PARAMETER (Z0ICE = 0.001)
PARAMETER (Z0SEA = 0.001)
*
*
*
*
#include "leads.cdk"
*
REAL ALMIN, TAUF, TAUDAY
*
*
DATA ALMIN / 0.50 /
DATA TAUF / 0.24 /
DATA TAUDAY / 24. /
*
SAVE ALMIN, TAUF, TAUDAY
*
*********************************************************************
*
*
*
EXTERNAL INISOILI
EXTERNAL CALCALB1, COHERENCE
*
************************************************************************
* AUTOMATIC ARRAYS
************************************************************************
*
AUTOMATIC ( ALDATD , REAL , (NCLASS))
AUTOMATIC ( D2DATD , REAL , (NCLASS))
AUTOMATIC ( RSMINDATD , REAL , (NCLASS))
AUTOMATIC ( LAIDATD , REAL , (NCLASS))
AUTOMATIC ( VEGDATD , REAL , (NCLASS))
AUTOMATIC ( CVDATD , REAL , (NCLASS))
AUTOMATIC ( RGLDATD , REAL , (NCLASS))
AUTOMATIC ( GAMMADATD , REAL , (NCLASS))
*
************************************************************************
*
#include "icelvls.cdk"
*
REAL prcor, diff
INTEGER I, K, M, IM
IM(I,M)=(M-1)*NI+I
*
****
*
*
*
*
* Several treatments on geophysical
* fields valid for both ISBA and FCREST
*
* The water temperature (TM) is decreased
* for points where the filtering of mountains
* lead to an icrease of the water level
* (old subroutine MODTMTP of GEM's dynamic library)
*
IF ( drylaps ) THEN
prcor = grav/cpd
ELSE
prcor = grav*stlo
END IF
*
*
*
*VDIR NODEP
DO i=0,ni-1
*
IF (montagn .AND. cortm) THEN
* MF and MT units are meters
diff = (e(mf+i) - e(mt+i)) * prcor
IF (cortm) THEN
IF (diff .GT. 0.) e(twateren+i) = e(twateren+i) - diff
END IF
END IF
*
*
IF (ischmsol .EQ. 1) THEN
IF (nint(e(veginden+i)) .EQ. 25) e(veginden+i) = 4.
IF (nint(e(veginden+i)) .EQ. 26) e(veginden+i) = 10.
END IF
*
END DO
*
*
* Other consistency tests ...
*
*VDIR NODEP
DO i=0,ni-1
e(snodpen +i ) = MAX( 0. , e(snodpen+i ))
END DO
*
*
*VDIR NODEP
DO i=0,ni-1
*
e(tglacen+i ) = MIN( TRPL, e(tglacen+i ))
e(tglacen+i+ni ) = MIN( TRPL, e(tglacen+i+ni ))
*
END DO
*
*
*
* From the "entry" to the "permanent" bus
*
*========================================================================
* for variables common to FCREST and ISBA
*========================================================================
*
*
*VDIR NODEP
DO i=0,ni-1
f(alvis +i+(indx_soil -1)*ni) = e(alen +i )
f(alvis +i+(indx_glacier-1)*ni) = e(alen +i )
f(alvis +i+(indx_water -1)*ni) = e(alen +i )
f(alvis +i+(indx_ice -1)*ni) = e(alen +i )
f(alvis +i+(indx_agrege -1)*ni) = e(alen +i )
f(dlat +i ) = e(dlaten +i )
f(dlon +i ) = e(dlonen +i )
f(glacier +i ) = e(glacen +i )
f(mg +i ) = e(mgen +i )
*
* --- snodp deja en metres
f(snodp +i+(indx_soil -1)*ni) = e(snodpen +i )
f(snodp +i+(indx_glacier-1)*ni) = e(snodpen +i )
f(snodp +i+(indx_water -1)*ni) = 0.0
f(snodp +i+(indx_ice -1)*ni) = e(snodpen +i )
*
f(twater +i ) = e(twateren +i )
f(tsoil +i ) = e(tsoilen +i )
f(tsoil +i + ni ) = e(tsoilen +i + ni)
f(tsrad +i ) = e(tsoilen +i )
f(z0 +i+(indx_soil -1)*ni) = e(z0en +i )
f(z0 +i+(indx_glacier-1)*ni) = z0ice
f(z0 +i+(indx_water -1)*ni) = z0sea
f(z0 +i+(indx_ice -1)*ni) = z0ice
f(z0 +i+(indx_agrege -1)*ni) = e(z0en +i )
f(z0t +i+(indx_soil -1)*ni) = e(z0en +i )
f(z0t +i+(indx_glacier-1)*ni) = z0ice
f(z0t +i+(indx_water -1)*ni) = z0sea
f(z0t +i+(indx_ice -1)*ni) = z0ice
f(z0t +i+(indx_agrege -1)*ni) = e(z0en +i )
f(lhtg +i ) = e(lhtgen +i )
f(icedp +i ) = e(icedpen +i )
f(tglacier+i ) = e(tglacen +i )
f(tglacier+i + ni ) = e(tglacen +i + ni)
f(glsea +i ) = e(glseaen +i )
f(glsea0 +i ) = e(glseaen +i )
* Mask for the lakes
f(ml +i ) = e(vegfen +i +2*ni)
*
* transvidage des variables necessaires au blocage orographique
f(dhdx +i) = e(dhdxen + i)
f(dhdy +i) = e(dhdyen + i)
f(dhdxdy +i) = e(dhdxdyen + i)
*
END DO
*
DO K=1,NL
DO I=0,NI-1
f( tmice +i + (k-1)*ni ) = e( tmicen +i + (k-1)*ni )
f( tmice +i + (k-1)*ni ) = min(tcdk, f( tmice +i + (k-1)*ni ))
END DO
END DO
*
* Special cases
*
call lacs(f, fsiz, ni)
*
*VDIR NODEP
do i=0,ni-1
* no snow allowed in the absence of marine ice
if (f(icedp+i).lt.himin) then
f(snodp+i+(indx_ice -1)*ni) = 0.0
endif
end do
*
*VDIR NODEP
DO i=0,ni-1
* For force-restore scheme only : if no radiation scheme
* is used, then surface IR emissivity is set to 0.
IF (iradia.EQ.0) THEN
f(epstfn+i) = 0.
ELSE
f(epstfn+i) = stefan
END IF
END DO
*
*
*========================================================================
* FOR FCREST ONLY
*========================================================================
*
IF ((schmsol.EQ.'FCREST').OR.(schmsol.EQ.'fcrest')) THEN
*
*
*VDIR NODEP
* Calculate the albedo for the FORCE-RESTORE
* scheme after saving original value in ALNOSNO.
DO i=0,ni-1
f(vegindx +i ) = e(veginden +i)
f(wsoil +i ) = e(hs +i)
f(wsoil +i +ni) = e(hs +i)
f(alnosno +i) = e(alen +i)
END DO
*
CALL calcalb1 (f(alnosno),f(vegindx),
+ f(snodp+(indx_soil-1)*ni),
+ f(alvis+(indx_soil-1)*ni), ni)
*
*
CALL coherence( f, fsiz, ni )
*
*
END IF
*
*=========================================================================
* FOR ISBA ... FOR ISBA ... FOR ISBA
*=========================================================================
*
*
IF ((schmsol.EQ.'ISBA').OR.(schmsol.EQ.'isba')) THEN
*
*
*
*VDIR NODEP
DO i=1,nclass*ni
*
f(vegf+i-1) = e(vegfen+i-1)
END DO
*
*VDIR NODEP
DO i=0,ni-1
*
f(wsoil +i ) = e(wsoilen +i )
f(wsoil +i+ni) = e(wsoilen +i + ni)
f(wveg +i ) = e(wvegen +i )
f(isoil +i ) = e(isoilen +i )
f(wsnow +i ) = e(wsnowen +i )
f(resa +i ) = 50.
END DO
*
*
* Special operations for the snow variables
*
* CAREFUL HERE about the units:
* "snoro" is the relative density of snow,
* i.e., rho_ice / rho_water (no units)
* "snoma" is the snow water equivalent in
* mm (i.e., kg / m2)
* "snoal" is the snow albedo determined from
* the snow age
*
* NOTE that "snoag" is in hours ... (tauday also)
*
*
*VDIR NODEP
DO i=0,ni-1
f(snoro + i) = MAX(100.,e(snoroen+i)) / rauw
f(snoma + i) = rauw * f(snoro+i) * f(snodp+i+(indx_soil-1)*ni)
END DO
*
*
* For the ALBEDO, there are two possibilities:
*
* 1) if switch "snoalb_anl" is true, then the "I6"
* record in the starting standard file (SNOALEN)
* contains the snow albedo
*
* 2) if switch "snoalb_anl" is false, then we
* use the snow age (SNOAGEN) to derive the
* snow albedo
*
IF (snoalb_anl) THEN
*
DO i=0,ni-1
f(snoal + i) = e(snoalen + i)
END DO
*
ELSE
*
* snow albedo is determined from the
* snow age according to two different
* expressions depending if the snow pack
* is melting or not
*
*VDIR NODEP
DO i=0,ni-1
IF (f(tmoins+ im(i,nk)).LT.0.) THEN
f(snoal + i) = ansmax - todry*e(snoagen+i)/tauday
ELSE
f(snoal + i) = (ansmax-almin) *
1 EXP( -tauf*e(snoagen+i)/tauday )
1 + almin
END IF
f(snoal + i) = MAX( f(snoal+i) , almin )
f(snoal + i) = MIN( f(snoal+i) , ansmax )
END DO
*
END IF
*
*
* Initialize the parameters that depend
* on vegetation
*
call iniveg( f, fsiz, 0, ni )
*
*
*
* Sand and clay fractions of the soil
* are taken as simple averages of the
* first 3 layers
*
*VDIR NODEP
DO i=0,ni-1
f(sand + i ) = ( e(sanden + i)
1 + e(sanden + i + ni)
1 + e(sanden + i + 2*ni) ) / 3.
f(clay + i ) = ( e(clayen + i)
1 + e(clayen + i + ni)
1 + e(clayen + i + 2*ni) ) / 3.
END DO
*
*
* Make sure the entry fields are
* coherent ...
*
CALL coherence( f, fsiz, ni )
*
*
* Initialize the soil characteristics
* using the soil texture
*
CALL inisoili( f, fsiz, ni )
*
*
END IF
*
*
*
RETURN
END