!copyright (C) 2001 MSC-RPN COMM %%%RPNPHY%%%
*** S/P SEAICE1
*
#include "phy_macros_f.h"
SUBROUTINE SEAICE1 ( BUS, BUSSIZ, 1,32
1 PTSURF, PTSURFSIZ,
1 TRNCH, KOUNT,
1 N, M, NK )
*
#include "impnone.cdk"
*
*
INTEGER BUSSIZ, KOUNT, TRNCH, N, M, NK
INTEGER PTSURFSIZ
INTEGER PTSURF(PTSURFSIZ)
REAL BUS(BUSSIZ)
*
*
*Author
* J. Mailhot (April 1999)
*
*Revisions
* 001 B. Bilodeau and S. Belair (April 1999)
* Coupling to the new SURFACE module
* 002 J. Mailhot (Feb 2000) - Add snow/ice melt/growth processes
* and revisions of some parameters. Change name of s/r.
* 003 B. Bilodeau (Nov 2000) - New comdeck sfcbus.cdk
* 004 B. Bilodeau (Jan 2001) - Automatic arrays
* 005 J. Mailhot (Apr 2001) - Multi-level ice model
* 006 B. Bilodeau (Aug 2001) - LOCBUS
* 007 J.-F. Mahfouf (Spring 2003) -
* Add implicit boundary condition option for vert. diff.
*Object
* The multi-level model calculates the temperature profile across a
* snow/ice slab and the surface fluxes of heat, moisture, and momentum
* over floating ice-covered surfaces (sea/lake).
*
*Arguments
*
* - Input/Output -
* BUS bus of surface variables
*
* - Input -
* BUSSIZ size of the surface bus
* PTSURF surface pointers
* PTSURFSIZ dimension of ptsurf
* KOUNT number of timestep
* DELT timestep
* ICEMELT switch to control ice and snow melting
* N running length
* M horizontal dimension
* NK vertical dimension
*
*
*Notes
* One-dimensional thermodynamic sea ice model:
* - based on modified version of nl-layer model of Semtner (1976, JPO 6, 379-389).
* - includes snow cover on the top of ice, heat conduction through snow and ice,
* thermal inertia of snow and ice layers, and penetrating solar radiation
* - includes a parameterization of albedo, conductivity and heat capacity
* (cf. Ebert and Curry 1993, JGR 98, 10085-10109; Flato and Brown 1996,
* JGR 101, 25767-25777)
* - effects of open leads in ice are considered (using climatological values
* of lead fraction to modify the analysed ice fraction)
* - vertical discretization similar to Flato and Brown (1996)
* but a flux boundary condition applied at the surface
*
* When the option ICEMELT = .TRUE. the model allows melting/growth of snow depth
* and ice thickness, and open water in summer with a crude oceanic mixed layer.
*
* The routine works with an arbitrary number of levels, but needs at least 2
* (NL .GE. 2).
**
*
************************************************************************
* AUTOMATIC ARRAYS
************************************************************************
*
#include "icelvls.cdk"
*
AUTOMATIC ( A , REAL , (N,NL) )
AUTOMATIC ( B , REAL , (N,NL) )
AUTOMATIC ( C , REAL , (N,NL) )
AUTOMATIC ( CAP , REAL , (N,NL) )
AUTOMATIC ( COND , REAL , (N,NL) )
AUTOMATIC ( D , REAL , (N,NL) )
AUTOMATIC ( DQSAT , REAL , (N ) )
AUTOMATIC ( DZ , REAL , (N,NL) )
AUTOMATIC ( EMIST , REAL , (N ) )
AUTOMATIC ( NSNOW , INTEGER , (N ) )
AUTOMATIC ( QSAT , REAL , (N ) )
AUTOMATIC ( RHOA , REAL , (N ) )
AUTOMATIC ( SCR1 , REAL , (N ) )
AUTOMATIC ( SCR2 , REAL , (N ) )
AUTOMATIC ( SCR3 , REAL , (N ) )
AUTOMATIC ( SCR4 , REAL , (N ) )
AUTOMATIC ( SCR5 , REAL , (N ) )
AUTOMATIC ( SCR6 , REAL , (N ) )
AUTOMATIC ( SCR7 , REAL , (N ) )
AUTOMATIC ( SCR8 , REAL , (N ) )
AUTOMATIC ( SCR9 , REAL , (N ) )
AUTOMATIC ( SCR10 , REAL , (N ) )
AUTOMATIC ( SCR11 , REAL , (N ) )
AUTOMATIC ( SOUR , REAL , (N,NL) )
AUTOMATIC ( TB , REAL , (N ) )
AUTOMATIC ( TP , REAL , (N,NL) )
AUTOMATIC ( VMOD , REAL , (N ) )
AUTOMATIC ( Z , REAL , (N,NL) )
AUTOMATIC ( ZTN , REAL , (N ) )
AUTOMATIC ( ZUN , REAL , (N ) )
*
************************************************************************
*
* Local variables "pointing" to the surface bus
*
REAL ALBSFC, CMU, CTU, FC_ICE
REAL FSOL, FV_ICE
REAL HICE, HST_ICE, HU, ILMO_ICE
REAL PS, QSICE, TH, TS, T, TT, UU, VV
REAL Z0H, Z0M
REAL ZALFAQ, ZALFAT, ZDLAT, ZFCOR, ZFDSI
REAL ZFTEMP, ZFVAP, ZQDIAG, ZSNODP, ZSNOWRATE, ZTDIAG
REAL ZTSRAD, ZUDIAG, ZVDIAG, ZUE2, ZZA
*
POINTER (IALBSFC , ALBSFC (1 ) )
POINTER (ICMU , CMU (1 ) )
POINTER (ICTU , CTU (1 ) )
POINTER (IFC_ICE , FC _ICE (1 ) )
POINTER (IFSOL , FSOL (1 ) )
POINTER (IFV_ICE , FV _ICE (1 ) )
POINTER (IHICE , HICE (1 ) )
POINTER (IHST_ICE , HST _ICE (1 ) )
POINTER (IHU , HU (1 ) )
POINTER (IILMO_ICE , ILMO _ICE (1 ) )
POINTER (IPS , PS (1 ) )
POINTER (IQSICE , QSICE (1 ) )
POINTER (IT , T (N,NL) )
POINTER (ITH , TH (1 ) )
POINTER (ITS , TS (1 ) )
POINTER (ITT , TT (1 ) )
POINTER (IUU , UU (1 ) )
POINTER (IVV , VV (1 ) )
POINTER (IZ0H , Z0H (1 ) )
POINTER (IZ0M , Z0M (1 ) )
POINTER (IZALFAQ , ZALFAQ (1 ) )
POINTER (IZALFAT , ZALFAT (1 ) )
POINTER (IZDLAT , ZDLAT (1 ) )
POINTER (IZFCOR , ZFCOR (1 ) )
POINTER (IZFDSI , ZFDSI (1 ) )
POINTER (IZFTEMP , ZFTEMP (1 ) )
POINTER (IZFVAP , ZFVAP (1 ) )
POINTER (IZQDIAG , ZQDIAG (1 ) )
POINTER (IZSNODP , ZSNODP (1 ) )
POINTER (IZSNOWRATE, ZSNOWRATE (1 ) )
POINTER (IZTDIAG , ZTDIAG (1 ) )
POINTER (IZTSRAD , ZTSRAD (1 ) )
POINTER (IZUDIAG , ZUDIAG (1 ) )
POINTER (IZVDIAG , ZVDIAG (1 ) )
POINTER (IZUE2 , ZUE2 (1 ) )
POINTER (IZZA , ZZA (1 ) )
*
*
INTEGER I, J, K
REAL BETA, SC
*
#include "vamin.cdk"
SAVE VAMIN
*
REAL CON1,CON2,CON3,CON4,CON5,CON6
REAL CON7,CON8,CON9,CON10,CON11,CON12
SAVE CON1,CON2,CON3,CON4,CON5,CON6
SAVE CON7,CON8,CON9,CON10,CON11,CON12
REAL FI0,CONDFI,TFRZW,TMELI,TMELS
SAVE FI0,CONDFI,TFRZW,TMELI,TMELS
REAL ALBOW,ALBDI,ALBMI,ALBDS,ALBMS,EMISI,EMISNO,EMISW
SAVE ALBOW,ALBDI,ALBMI,ALBDS,ALBMS,EMISI,EMISNO,EMISW
REAL COEFCOND,COEFHCAP,COEFEXT
SAVE COEFCOND,COEFHCAP,COEFEXT
REAL ROICE,ROSNOW(2),ROSWTR
SAVE ROICE,ROSNOW,ROSWTR
REAL Z0ICE,Z0W
SAVE Z0ICE,Z0W
REAL BASEHF
SAVE BASEHF
REAL HCAPI,HCAPW,VHFICE,VHFBAS,VHFSNO
SAVE HCAPI,HCAPW,VHFICE,VHFBAS,VHFSNO
REAL HSMIN,DMIX
SAVE HSMIN,DMIX
*
DATA CON1 , CON2 , CON3 , CON4 /
1 2.845E-6 , 2.7E-4 , 233.0 , 0.2 /
DATA CON5 , CON6 , CON7 , CON8 /
1 92.88 , 7.364 , 3.2 , 14.24 /
DATA CON9 , CON10 , CON11 , CON12 /
1 19.39 , 0.1 , 0.44 , 0.075 /
DATA TFRZW , TMELI , TMELS / 271.2 , 273.05 , 273.15 /
DATA ALBOW , ALBDI , ALBMI , ALBDS , ALBMS /
1 0.08 , 0.57 , 0.50 , 0.83 , 0.77 /
DATA FI0 , CONDFI , COEFCOND , COEFHCAP , COEFEXT /
1 0.17 , 2.034 , 0.1172 , 1.715E+7 , 1.5 /
DATA EMISI , EMISNO , EMISW / 0.99 , 0.99 , 0.97 /
DATA ROICE , ROSWTR / 913.0 , 1025.0 /
DATA ROSNOW / 330.0 , 450.0 /
DATA Z0ICE , Z0W / 1.6E-4 , 3.2E-5 /
DATA BASEHF / 2.0 /
DATA HCAPI , HCAPW , VHFICE , VHFBAS , VHFSNO /
1 2.062E+3 , 4.088E+3 , 3.014E+8 , 2.679E+8 , 1.097E+8 /
DATA HSMIN , DMIX / 0.010 , 30.0 /
*
#include "himin.cdk"
*
INTEGER PTR, X
*
#include "locbus.cdk"
INTEGER INDX_SFC, SURFLEN
PARAMETER (INDX_SFC = INDX_ICE)
INTEGER QUELNIVO(MAXVARSURF)
*
#include "zuzt.cdk"
*
#include "consphy.cdk"
*
#include "options.cdk"
*
#include "dintern.cdk"
*
#include "sfcbus.cdk"
*
* fonctions-formule
*
#include "xptsurf.cdk"
*
#include "fintern.cdk"
*
*
*
SURFLEN = M
*
*
* EQUIVALENCES
*
INIT_LOCBUS()
*
* Syntax of macro locbus (must be typed in CAPITAL letters):
* locbus (pointer, array_name_in_the_bus, level)
* If level=0, array chosen automatically as follows:
* 1) level = 1 if array has 1 level only (e.g. TSURF )
* 2) level = nk if array has nk levels (e.g. TMOINS)
* 3) level = indx_sfc if array has a level for each surface type (e.g. FC)
* 4) level has to be specified by user if array has more than one level
* that all "belong" to the same surface type (e.g. TSOIL)
*
LOCBUS (IALBSFC , ALVIS , 0 )
LOCBUS (ICMU , BM , 0 )
LOCBUS (ICTU , BT , 0 )
LOCBUS (IFC _ICE , FC , 0 )
LOCBUS (IFSOL , FLUSOLIS, 0 )
LOCBUS (IFV _ICE , FV , 0 )
LOCBUS (IHICE , ICEDP , 0 )
LOCBUS (IHST _ICE , HST , 0 )
LOCBUS (IHU , HUMOINS , 0 )
LOCBUS (IILMO_ICE , ILMO , 0 )
LOCBUS (IPS , PMOINS , 0 )
LOCBUS (IQSICE , QSURF , 0 )
LOCBUS (IT , TMICE , 1 )
LOCBUS (ITH , THETAA , 0 )
LOCBUS (ITS , TSURF , 0 )
LOCBUS (ITT , TMOINS , 0 )
LOCBUS (IUU , UMOINS , 0 )
LOCBUS (IVV , VMOINS , 0 )
LOCBUS (IZ0H , Z0T , 0 )
LOCBUS (IZ0M , Z0 , 0 )
LOCBUS (IZALFAQ , ALFAQ , 0 )
LOCBUS (IZALFAT , ALFAT , 0 )
LOCBUS (IZDLAT , DLAT , 0 )
LOCBUS (IZFCOR , FCOR , 0 )
LOCBUS (IZFDSI , FDSI , 0 )
LOCBUS (IZFTEMP , FTEMP , 0 )
LOCBUS (IZFVAP , FVAP , 0 )
LOCBUS (IZSNODP , SNODP , 0 )
LOCBUS (IZSNOWRATE, SNOWRATE, 0 )
LOCBUS (IZTSRAD , TSRAD , 0 )
LOCBUS (IZUDIAG , UDIAG , 0 )
LOCBUS (IZVDIAG , VDIAG , 0 )
LOCBUS (IZTDIAG , TDIAG , 0 )
LOCBUS (IZQDIAG , QDIAG , 0 )
LOCBUS (IZUE2 , UE2 , 0 )
LOCBUS (IZZA , ZA , 0 )
*
*
DO I=1,N
ZTN(I) = ZT
ZUN(I) = ZU
END DO
*
* test on minimum value for number of layers
IF( NL .LT. 2 ) THEN
PRINT *,'******* FATAL ERROR IN ICE MODULE - NL < 2 *******'
STOP
ENDIF
*
* fully-implicit time scheme
BETA = 1.0
SC = (1.0 - BETA)*DELT
*
*
** 1. Preliminaries
* --------------------
*
DO I=1,N
*
* Air density near the surface
RHOA(I) = PS(I)/(RGASD*TT(I)*(1.0+DELTA*HU(I)))
*
* Modifications for the sea ice surface:
* - under-ice seawater temperature
TB(I) = TFRZW
*
IF( ICEMELT ) THEN
* Roughness lengths for the surface (ice/water)
IF( HICE(I) .GE. HIMIN ) THEN
Z0M(I) = Z0ICE
ELSE
Z0M(I) = Z0W
* remove snow if ice is too thin
ZSNODP(I) = 0.0
ENDIF
ELSE
* - minimum ice thickness
HICE(I) = MAX ( HICE(I) , HIMIN )
Z0M(I) = Z0ICE
*
ENDIF
*
Z0H(I) = Z0M(I)
* Wind module
*
VMOD (I) = SQRT( MAX( VAMIN,(UU(I)**2 + VV(I)**2)))
*
END DO
*
*
* 2. Initialization of temperature profiles
* ---------------------------------------------
*
*
* Temperature at mid-layers
DO K=1,NL-1
DO I=1,N
TP(I,K) = 0.5*( T(I,K) + T(I,K+1) )
END DO
END DO
*
DO I=1,N
TP(I,NL) = 0.5*( T(I,NL) + TB(I) )
END DO
*
*
* 3. Calculate the drag and heat coefficients
* -----------------------------------------------
*
* Arbitrary values for the height
* of the boundary layer (for FLXSURF1)
* No impact on the results.
DO I=1,N
* Saturated specific humidity at surface
TS(I) = T(I,1)
QSAT(I) = FOQST ( TS (I), PS(I) )
DQSAT(I) = FODQS ( QSAT(I), TS(I) )
END DO
*
CALL FLXSURF2( CMU, CTU, SCR1, ZFTEMP, ZFVAP, ILMO_ICE,
$ ZUE2, ZFCOR, TH, HU, ZZA, VMOD, TS,
$ QSICE, HST_ICE, Z0M, Z0H,
$ SCR2, SCR3, SCR4, SCR5, SCR6, SCR7, N)
*
*
*
* 4. Parameterizations (albedo, conductivity, capacity,...)
* -------------------------------------------------------------
*
*
DO I=1,N
* - surface albedo (function of surface type and temperature)
IF( TS(I) .LT. TMELS ) THEN
SCR2(I) = ALBDS
ELSE
SCR2(I) = ALBMS
ENDIF
*
IF( TS(I) .LT. TMELI ) THEN
SCR3(I) = MIN(ALBDI,MAX(ALBOW,0.08+CON11*HICE(I)**0.28))
SCR9(I) = CHLC+CHLF
ELSE
SCR3(I) = ALBMI
SCR9(I) = CHLC
ENDIF
*
END DO
*
DO I=1,N
* - emissivity and fraction of penetrating solar radiation
* (function of surface type snow/ice/water)
IF( HICE(I) .GE. HIMIN ) THEN
IF( ZSNODP(I) .GT. CON10 ) THEN
ALBSFC(I) = SCR2(I)
EMIST(I) = EMISNO
SCR6(I) = 1.0
ELSE
ALBSFC(I) = MIN( SCR2(I) , SCR3(I)+ZSNODP(I)*
1 (SCR2(I)-SCR3(I))/CON10 )
EMIST(I) = EMISI
SCR6(I) = 1.0 - FI0
SCR6(I) = MIN( 1.0 , SCR6(I)+ZSNODP(I)*
1 (1.0-SCR6(I))/CON10 )
ENDIF
ELSE
ALBSFC(I) = ALBOW
EMIST(I) = EMISW
SCR6(I) = 1.0
SCR9(I) = CHLC
ENDIF
* - bulk salinity of the ice
SCR8(I) = MAX( CON7 , CON8-CON9*HICE(I) )
*
END DO
*
*
DO I=1,N
IF( HICE(I) .GE. HIMIN ) THEN
* snow layer
IF( ZSNODP(I) .GE. HSMIN ) THEN
NSNOW(I) = 1
DZ(I,1) = ZSNODP(I)
SCR4(I) = ROSNOW(1)
IF( TP(I,1) .GE. TMELS ) SCR4(I) = ROSNOW(2)
COND(I,1) = CON1*SCR4(I)**2+CON2*2.**((TP(I,1)-CON3)*CON4)
CAP(I,1) = SCR4(I)*(CON5+CON6*TP(I,1))
ELSE
NSNOW(I) = 0
DZ(I,1) = HICE(I)/FLOAT(NL)
ENDIF
*
ELSE
NSNOW(I) = 0
DZ(I,1) = DMIX
CAP(I,1) = DMIX*ROSWTR*HCAPW
ENDIF
END DO
*
* ice slab
DO K=2,NL
DO I=1,N
IF( HICE(I) .GE. HIMIN ) THEN
SCR3(I) = ROICE*HCAPI
DZ(I,K) = HICE(I)/FLOAT(NL-NSNOW(I))
COND(I,K) = CONDFI+COEFCOND*SCR8(I)/MIN(TP(I,K)-TMELI,-0.1)
COND(I,K) = MAX( 0.2*CONDFI , COND(I,K) )
*
CAP(I,K) = SCR3(I)+COEFHCAP*SCR8(I)/
1 MIN( T(I,K)-TMELI , -0.1 )**2
CAP(I,K) = MIN( 100.0*SCR3(I) , CAP(I,K) )
ELSE
NSNOW(I) = 0
DZ(I,K) = DMIX
CAP(I,K) = DMIX*ROSWTR*HCAPW
ENDIF
END DO
END DO
*
DO I=1,N
IF( HICE(I) .GE. HIMIN .AND. NSNOW(I).EQ.0 ) THEN
SCR3(I) = ROICE*HCAPI
DZ(I,1) = HICE(I)/FLOAT(NL-NSNOW(I))
COND(I,1) = CONDFI+COEFCOND*SCR8(I)/MIN(TP(I,1)-TMELI,-0.1)
COND(I,1) = MAX( 0.2*CONDFI , COND(I,1) )
*
CAP(I,1) = SCR3(I)+COEFHCAP*SCR8(I)/
1 MIN( T(I,1)-TMELI , -0.1 )**2
CAP(I,1) = MIN( 100.0*SCR3(I) , CAP(I,1) )
ENDIF
END DO
*
DO I=1,N
IF( HICE(I) .GE. HIMIN ) THEN
IF( NSNOW(I) .EQ. 0 ) THEN
CAP(I,1) = SCR3(I)+COEFHCAP*SCR8(I)/
1 MIN( 0.75*T(I,1)+0.25*T(I,2)-TMELI , -0.1 )**2
CAP(I,1) = MIN( 100.0*SCR3(I) , CAP(I,1) )
ELSE
IF( NL .GT. 2) THEN
CAP(I,2) = SCR3(I)+COEFHCAP*SCR8(I)/
1 MIN( 0.75*T(I,2)+0.25*T(I,3)-TMELI , -0.1 )**2
ELSE
CAP(I,2) = SCR3(I)+COEFHCAP*SCR8(I)/
1 MIN( 0.75*T(I,2)+0.25*TB(I)-TMELI , -0.1 )**2
ENDIF
CAP(I,2) = MIN( 100.0*SCR3(I) , CAP(I,2) )
ENDIF
* add "thin" snow layer
IF( NSNOW(I) .EQ. 0 ) THEN
SCR4(I) = ROSNOW(1)
IF( T(I,1) .GE. TMELS ) SCR4(I) = ROSNOW(2)
SCR4(I) = CON1*SCR4(I)**2+CON2*2.**((T(I,1)-CON3)*CON4)
COND(I,1) = COND(I,1)/
1 (1.0+ZSNODP(I)*COND(I,1)/(DZ(I,1)*SCR4(I)))
ENDIF
*
ENDIF
END DO
*
DO I=1,N
* penetrating solar radiation
SCR7(I) = (1.0-SCR6(I))*FSOL(I)*(1.-ALBSFC(I))
IF( NSNOW(I).EQ.1 ) THEN
Z(I,1) = 0.0
Z(I,2) = 0.5*DZ(I,2)
ELSE
Z(I,1) = 0.5*DZ(I,1)
Z(I,2) = Z(I,1)+0.5*(DZ(I,2)+DZ(I,1))
ENDIF
SOUR(I,1) = SCR7(I)*(1.-EXP(-COEFEXT*Z(I,1)))
SOUR(I,2) = SCR7(I)*(EXP(-COEFEXT*Z(I,1))
1 -EXP(-COEFEXT*Z(I,2)))
*
END DO
*
IF( NL .GT. 2) THEN
DO K=3,NL
DO I=1,N
Z(I,K) = Z(I,K-1)+0.5*(DZ(I,K)+DZ(I,K-1))
SOUR(I,K) = SCR7(I)*(EXP(-COEFEXT*Z(I,K-1))
1 -EXP(-COEFEXT*Z(I,K)))
END DO
END DO
DO I=1,N
Z(I,NL) = Z(I,NL)+0.5*DZ(I,NL)
SOUR(I,NL) = SCR7(I)*(EXP(-COEFEXT*Z(I,NL-1))
1 -EXP(-COEFEXT*Z(I,NL)))
END DO
ELSE
DO I=1,N
Z(I,NL) = Z(I,NL)+0.5*DZ(I,NL)
SOUR(I,NL) = SCR7(I)*(1.0-EXP(-COEFEXT*Z(I,NL)))
END DO
ENDIF
*
*
*
* 5. Compute the temperature profile
* --------------------------------------
*
* linearized terms in surface heat budget
DO I=1,N
*
SCR1(I) = 4. * EMIST(I) * STEFAN * TS(I)**3
1 + RHOA(I) * CTU(I) * (DQSAT(I) * SCR9(I) + CPD)
1
*
SCR2(I) = 3. * EMIST(I) * STEFAN * TS(I)**4
1 + RHOA(I) * CTU(I) * DQSAT(I) * SCR9(I) * TS(I)
*
SCR3(I) = RHOA(I)*CTU(I)*(CPD*TH(I)-SCR9(I)*(QSAT(I)-HU(I)))
1 + SCR6(I)*FSOL(I)*(1.-ALBSFC(I))
1 + EMIST(I)*ZFDSI(I)
*
END DO
*
* setup tridiagonal terms A, B, C
* and right-hand-side term D
*
IF( NL.GT.3 ) THEN
DO K=3,NL-1
DO I=1,N
IF( HICE(I) .GE. HIMIN )
1 CAP(I,K) = 0.5*( DZ(I,K)+DZ(I,K-1) )*CAP(I,K)
END DO
END DO
ENDIF
*
DO I=1,N
IF( HICE(I) .GE. HIMIN ) THEN
SOUR(I,2) = SOUR(I,2)+SOUR(I,1)
IF( NSNOW(I) .EQ. 0 ) THEN
CAP(I,2) = ( 0.5*DZ(I,2)+DZ(I,1) )*CAP(I,2)
* add "thin" snow layer
SCR4(I) = ROSNOW(1)
IF( T(I,1) .GE. TMELS ) SCR4(I) = ROSNOW(2)
SCR4(I) = SCR4(I)*(CON5+CON6*SCR5(I))
CAP(I,2) = CAP(I,2) + ZSNODP(I)*SCR4(I)
CAP(I,NL) = 0.5*( DZ(I,NL-1)+DZ(I,NL) )*CAP(I,NL)
ELSE
CAP(I,2) = 0.5*DZ(I,2)*CAP(I,2)+DZ(I,1)*CAP(I,1)
CAP(I,NL) = 0.5*( DZ(I,NL-1)+DZ(I,NL) )*CAP(I,NL)
ENDIF
ENDIF
END DO
*
DO K=2,NL
DO I=1,N
IF( HICE(I) .GE. HIMIN ) THEN
A(I,K) = (COND(I,K-1)/DZ(I,K-1))/CAP(I,K)
C(I,K) = (COND(I,K)/DZ(I,K))/CAP(I,K)
B(I,K) = -A(I,K)-C(I,K)
D(I,K) = DELT*SOUR(I,K)/CAP(I,K) + T(I,K)
ENDIF
END DO
END DO
*
DO I=1,N
IF( HICE(I) .GE. HIMIN ) THEN
A(I,1) = 0.0
B(I,1) = 0.0
C(I,1) = 0.0
D(I,1) = 0.0
C(I,NL) = 0.0
ENDIF
END DO
*
DO K=1,NL
DO I=1,N
IF( HICE(I) .LT. HIMIN ) THEN
A(I,K) = 0.0
B(I,K) = 0.0
C(I,K) = 0.0
D(I,K) = 0.0
ENDIF
END DO
END DO
*
DO I=1,N
IF( HICE(I) .LT. HIMIN )
1 D(I,1) = T(I,1)
END DO
*
*
*
* back to full levels
DO K=1,NL
DO I=1,N
TP(I,K) = T(I,K)
END DO
END DO
*
CALL DIFUVD1(D, SC, A, B, C, TP, D, N, N, NL)
*
DO K=1,NL
DO I=1,N
A(I,K) = -BETA*DELT*A(I,K)
B(I,K) = 1.0-BETA*DELT*B(I,K)
C(I,K) = -BETA*DELT*C(I,K)
END DO
END DO
*
DO I=1,N
* add upper boundary condition
IF( HICE(I) .GE. HIMIN ) THEN
C(I,1) = -BETA
B(I,1) = -C(I,1)+DZ(I,1)*SCR1(I)/COND(I,1)
D(I,1) = DZ(I,1)*(SCR2(I)+SCR3(I))/COND(I,1)
1 +(1.0-BETA)*(T(I,2)-T(I,1))
* add lower boundary condition
D(I,NL) = D(I,NL)+DELT*TB(I)*(COND(I,NL)/
1 DZ(I,NL))/CAP(I,NL)
ENDIF
END DO
*
DO K=1,NL
DO I=1,N
IF( HICE(I) .LT. HIMIN ) THEN
A(I,K) = -1.0
B(I,K) = 1.0
ENDIF
END DO
END DO
*
DO I=1,N
IF( HICE(I) .LT. HIMIN ) THEN
A(I,1) = 0.0
B(I,1) = 1.0
C(I,1) = 0.0
B(I,1) = B(I,1)+DELT*SCR1(I)/CAP(I,1)
D(I,1) = D(I,1)+DELT*(SCR2(I)+SCR3(I))/CAP(I,1)
ENDIF
END DO
*
*
CALL DIFUVD2(TP, A, B, C, D, D, N, N, NL)
*
* prevent temperatures
* from exceeding melting temperature
DO I=1,N
IF( HICE(I) .GE. HIMIN ) THEN
IF( ZSNODP(I) .GE. HSMIN ) THEN
SCR4(I) = TMELS
ELSE
SCR4(I) = TMELI
ENDIF
ENDIF
END DO
*
DO K=1,NL
DO I=1,N
IF( HICE(I) .GE. HIMIN )
1 TP(I,K) = MIN ( TP(I,K) , SCR4(I) )
END DO
END DO
*
DO I=1,N
IF( HICE(I) .LT. HIMIN )
1 SCR4(I) = TB(I)
END DO
*
*
* surface and snow/ice temperatures
DO K=1,NL
DO I=1,N
T(I,K) = TP(I,K)
END DO
END DO
*
DO I=1,N
TS(I) = TP(I,1)
* saturated specific humidity at surface
QSICE(I) = FOQST( TS(I), PS(I) )
END DO
*
*
* 6. Melting and growth
* -------------------------
* melting at the surface (snow/ice)
* growth/melting of ice at the lower boundary
IF(ICEMELT) THEN
*
DO I=1,N
*
* recompute surface albedo with new TS
*
IF( TS(I) .LT. TMELS ) THEN
SCR2(I) = ALBDS
ELSE
SCR2(I) = ALBMS
ENDIF
*
IF( TS(I) .LT. TMELI ) THEN
SCR3(I) = MIN(ALBDI,MAX(ALBOW,0.08+CON11*HICE(I)**0.28))
SCR9(I) = CHLC+CHLF
ELSE
SCR3(I) = ALBMI
SCR9(I) = CHLC
ENDIF
*
IF( HICE(I) .GE. HIMIN ) THEN
IF( ZSNODP(I) .GT. CON10 ) THEN
ALBSFC(I) = SCR2(I)
ELSE
ALBSFC(I) = MIN( SCR2(I) , SCR3(I)+ZSNODP(I)*
1 (SCR2(I)-SCR3(I))/CON10 )
ENDIF
ELSE
ALBSFC(I) = ALBOW
SCR9(I) = CHLC
ENDIF
*
END DO
*
DO I=1,N
* compute heat conduction flux in upper layer
IF( HICE(I) .GE. HIMIN ) THEN
IF( ZSNODP(I) .GE. HSMIN ) THEN
SCR7(I) = ROSNOW(1)
SCR5(I) = 0.5*(T(I,1)+T(I,2))
IF( SCR5(I) .GE. TMELS ) SCR7(I) = ROSNOW(2)
COND(I,1) = CON1*SCR7(I)**2
1 +CON2*2.**((SCR5(I)-CON3)*CON4)
ELSE
SCR7(I) = 0.5*(T(I,1)+T(I,2))
COND(I,1) = CONDFI+COEFCOND*SCR8(I)
1 /MIN(SCR7(I)-TMELI,-0.1)
COND(I,1) = MAX( 0.2*CONDFI , COND(I,1) )
* add "thin" snow layer
SCR5(I) = ROSNOW(1)
IF( T(I,1) .GE. TMELS ) SCR5(I) = ROSNOW(2)
SCR5(I) = CON1*SCR5(I)**2+CON2*2.**((T(I,1)-CON3)*CON4)
COND(I,1) = COND(I,1)/
1 (1.0+ZSNODP(I)*COND(I,1)/(DZ(I,1)*SCR5(I)))
ENDIF
SCR11(I) = COND(I,1)*(T(I,2)-T(I,1))/DZ(I,1)
ELSE
SCR11(I) = 0.0
ENDIF
*
* compute HA = H + LE + FWnet + FLnet
SCR5(I) = RHOA(I)*CTU(I)*
1 ( CPD*(TS(I)-TH(I))+SCR9(I)*(QSICE(I)-HU(I)) )
1 - SCR6(I)*FSOL(I)*(1.-ALBSFC(I))
1 + EMIST(I)*( STEFAN*TS(I)**4 - ZFDSI(I) )
*
SCR3(I) = SCR5(I) - SCR11(I)
*
END DO
*
*
DO I=1,N
* criteria for melting at the surface
SCR10(I) = 0.0
IF( TS(I).GE.SCR4(I) .AND.
1 SCR5(I).LE.0.0 .AND. SCR3(I).LT.0.0 ) THEN
* melt available snow...
ZSNODP(I) = ZSNODP(I) + DELT*SCR3(I)/VHFSNO
* ...then ice...
IF( ZSNODP(I).LT.0.0 ) THEN
HICE(I) = HICE(I) + ZSNODP(I)*VHFSNO/VHFICE
ZSNODP(I) = 0.0
* ...then warm oceanic mixed layer
IF( HICE(I).LT.0.0 ) THEN
TS(I) = TS(I) - HICE(I)*VHFICE/CAP(I,1)
HICE(I) = 0.0
SCR10(I) = 1.0
ENDIF
ENDIF
*
ENDIF
END DO
* add growth/melt of ice at lower boundary
DO I=1,N
IF( HICE(I) .GE. HIMIN ) THEN
*
* compute heat conduction flux in lower layer
SCR11(I) = COND(I,NL)*(TB(I)-T(I,NL))/DZ(I,NL)
* compute absorption of penetrating solar radiation
SCR8(I) = FSOL(I)*(1.-ALBSFC(I))*(1.-SCR6(I))
1 *EXP(-COEFEXT*Z(I,NL))
HICE(I) = MAX ( 0.0 , HICE(I) +
1 DELT*( SCR11(I)-SCR8(I)-BASEHF )/VHFBAS )
* if ice is too thin, switch to oceanic mixed layer
IF( HICE(I).LT.HIMIN ) THEN
TS(I) = TFRZW
ZSNODP(I) = 0.0
SCR10 (I) = 1.0
ENDIF
*
ELSE
IF( SCR5(I).GT.0.0 ) THEN
* cool oceanic mixed layer...
TS(I) = TS(I) - DELT*SCR5(I)/CAP(I,1)
* ...then form new ice
IF( TS(I).LE.TFRZW ) THEN
HICE(I) = HICE(I) + (TFRZW-TS(I))*CAP(I,1)/VHFBAS
TS(I) = TFRZW
ENDIF
SCR10(I) = 1.0
ENDIF
*
ENDIF
END DO
* ...and adjust final temperature profile
DO K=1,NL
DO I=1,N
IF( SCR10(I).EQ.1.0 ) THEN
T(I,K) = TS(I)
TP(I,K) = TS(I)
ENDIF
END DO
END DO
* add snow fall (units changed from
* water equivalent to snow equivalent)
DO I=1,N
IF( TS(I).LT.TMELI .AND. HICE(I).GT.HIMIN ) THEN
ZSNODP(I) = ZSNODP(I) + DELT*ZSNOWRATE(I)*(1000./ROSNOW(1))
ENDIF
END DO
*
* transition periods
DO I=1,N
IF( NSNOW(I).EQ.0 .AND. ZSNODP(I).GE.HSMIN ) THEN
T(I,1) = TP(I,1)
T(I,2) = TP(I,1)
Z(I,1) = 0.0
ENDIF
END DO
DO K=2,NL
DO I=1,N
IF( NSNOW(I).EQ.0 .AND. ZSNODP(I).GE.HSMIN )
1 Z(I,K) = Z(I,K-1)+DZ(I,K-1)
END DO
END DO
*
DO I=1,N
IF( NSNOW(I).EQ.0 .AND. ZSNODP(I).GE.HSMIN ) THEN
DZ(I,1) = 0.0
DZ(I,2) = 0.0
ENDIF
END DO
*
IF( NL.GT.2 ) THEN
DO K=3,NL
DO I=1,N
IF( NSNOW(I).EQ.0 .AND. ZSNODP(I).GE.HSMIN ) THEN
DZ(I,K) = DZ(I,K-1)+HICE(I)/FLOAT(NL-1)
T(I,K) = TP(I,K-1)+((DZ(I,K)-Z(I,K-1))
1 /(Z(I,K)-Z(I,K-1)))*(TP(I,K)-TP(I,K-1))
ENDIF
END DO
END DO
ENDIF
*
DO I=1,N
IF( NSNOW(I).EQ.1 .AND. ZSNODP(I).LT.HSMIN ) THEN
T(I,1) = TP(I,2)
Z(I,1) = 0.0
Z(I,2) = 0.0
ENDIF
END DO
IF( NL.GT.2 ) THEN
DO K=3,NL
DO I=1,N
IF( NSNOW(I).EQ.1 .AND. ZSNODP(I).LT.HSMIN )
1 Z(I,K) = Z(I,K-1)+DZ(I,K-1)
END DO
END DO
ENDIF
DO I=1,N
IF( NSNOW(I).EQ.1 .AND. ZSNODP(I).LT.HSMIN ) THEN
DZ(I,1) = 0.0
DZ(I,2) = HICE(I)/FLOAT(NL)
T(I,2) = TP(I,2)+((DZ(I,2)-Z(I,2))
1 /(Z(I,3)-Z(I,2)))*(TP(I,3)-TP(I,2))
ENDIF
END DO
IF( NL.GT.3 ) THEN
DO K=3,NL-1
DO I=1,N
IF( NSNOW(I).EQ.1 .AND. ZSNODP(I).LT.HSMIN ) THEN
DZ(I,K) = DZ(I,K-1)+HICE(I)/FLOAT(NL)
T(I,K) = TP(I,K)+((DZ(I,K)-Z(I,K))
1 /(Z(I,K+1)-Z(I,K)))*(TP(I,K+1)-TP(I,K))
ENDIF
END DO
END DO
ENDIF
DO I=1,N
IF( NSNOW(I).EQ.1 .AND. ZSNODP(I).LT.HSMIN ) THEN
DZ(I,NL) = DZ(I,NL-1)+HICE(I)/FLOAT(NL)
T(I,NL) = TP(I,NL)+((DZ(I,NL)-Z(I,NL))
1 /(MAX(HICE(I),HIMIN)/FLOAT(NL-1)))*(TB(I)-TP(I,NL))
ENDIF
END DO
*
ENDIF
*
*
* 7. Update fluxes
* ------------------------------------------------
*
* Update TS and QSICE
DO I=1,N
TS (I) = T(I,1)
QSICE(I) = FOQST( TS(I), PS(I) )
END DO
*
CALL DIASURF1(ZUDIAG, ZVDIAG, ZTDIAG, ZQDIAG,
$ N, UU, VV, TS, QSICE,
$ Z0M, Z0H, ILMO_ICE, ZZA,
$ HST_ICE, ZUE2, ZFTEMP, ZFVAP,
$ ZUN, ZTN, ZDLAT)
*
*
*VDIR NODEP
DO I=1,N
* remove snow if ice is too thin
IF( HICE(I).LT.HIMIN ) ZSNODP(I) = 0.0
*
ZTSRAD (I) = TS(I)
*
ZALFAT (I) = - CTU(I) * ( TS (I) - TH(I) )
ZALFAQ (I) = - CTU(I) * ( QSICE (I) - HU(I) )
IF (.NOT.IMPFLX) CTU (I) = 0.
RHOA (I) = PS(I)/(RGASD * ZTDIAG(I)*(1.+DELTA*ZQDIAG(I)))
FC _ICE(I) = -CPD *RHOA(I)*ZALFAT(I)
FV _ICE(I) = -(CHLC+CHLF)*RHOA(I)*ZALFAQ(I)
IF (IMPFLX) THEN
ZALFAT (I) = - CTU(I) * TS(I)
ZALFAQ (I) = - CTU(I) * QSICE(I)
ENDIF
****
END DO
*
* FILL THE ARRAYS TO BE AGGREGATED LATER IN S/R AGREGE
CALL FILLAGG ( BUS, BUSSIZ, PTSURF, PTSURFSIZ, INDX_ICE,
+ SURFLEN )
*
RETURN
END