MODULE STATIC_OCEAN 9,9
!@sum STATIC_OCEAN contains the ocean subroutines common to all Q-flux
!@+ and fixed SST runs
!@auth Original Development Team
!@ver 1.0 (Q-flux ocean)
!@cont OSTRUC,OCLIM,init_OCEAN,daily_OCEAN,DIAGCO
!@+ PRECIP_OC,OCEANS
USE CONSTANT
, only : lhm,rhows,rhoi,shw,shi,by12,byshi
USE MODEL_COM, only : im,jm,lm,focean,fland,fearth,flice
* ,Iyear1,Itime,jmon,jdate,jday,jyear,jmpery,JDendOfM,JDmidOfM
* ,ItimeI,kocean,itocean,itoice
USE GEOM
USE PBLCOM, only : npbl,uabl,vabl,tabl,qabl,eabl,cm=>cmgs,ch=>chgs
* ,cq=>cqgs,ipbl,roughl
USE SEAICE, only : xsi,ace1i,z1i,ac2oim,z2oim,ssi0,tfrez,fleadoc
* ,lmi
USE SEAICE_COM, only : rsi,msi,hsi,snowi,ssi
USE LANDICE_COM, only : snowli,tlandi
USE FLUXES, only : gtemp,sss,fwsim,mlhc
USE DAGCOM, only : aij, areg, jreg, ij_smfx, aj, j_implh, j_implm,
* j_imelt, j_hmelt, j_smelt
IMPLICIT NONE
SAVE
logical :: off_line=.false.
!@dbparam ocn_cycl determines whether prescribed ocean data repeat
!@+ after 1 year - 0:no 1:yes 2:seaice repeats,sst does not
integer :: ocn_cycl = 1
!@var TOCEAN temperature of the ocean (C)
REAL*8, DIMENSION(3,IM,JM) :: TOCEAN
!@var SSS0 default sea surface salinity (psu)
REAL*8 :: SSS0=34.7d0
!@var OTA,OTB,OTC ocean heat transport coefficients
REAL*8, DIMENSION(IM,JM,4) :: OTA,OTB
REAL*8, DIMENSION(IM,JM) :: OTC
!@var SINANG,COSANG Fourier coefficients for heat transports
REAL*8 :: SINANG,SN2ANG,SN3ANG,SN4ANG,
* COSANG,CS2ANG,CS3ANG,CS4ANG
!@var Z1O ocean mixed layer depth
REAL*8, DIMENSION(IM,JM) :: Z1O
!@var Z1OOLD previous ocean mixed layer depth
REAL*8, DIMENSION(IM,JM) :: Z1OOLD
!@var Z12O annual maximum ocean mixed layer depth
REAL*8, DIMENSION(IM,JM) :: Z12O
REAL*8, DIMENSION(IM,JM) :: DM,AOST,EOST1,EOST0,BOST
* ,COST,ARSI,ERSI1,ERSI0,BRSI,CRSI
INTEGER, DIMENSION(IM,JM) :: KRSI
COMMON/OOBS/DM,AOST,EOST1,EOST0,BOST,COST,ARSI,ERSI1,ERSI0,BRSI
* ,CRSI,KRSI
!@var iu_OSST,iu_SICE,iu_OCNML unit numbers for climatologies
INTEGER iu_OSST,iu_SICE,iu_OCNML
CONTAINS
SUBROUTINE OSTRUC(QTCHNG) 1
!@sum OSTRUC restructures the ocean temperature profile as ML
!@sum depths are changed (generally once a day)
!@auth Original Development Team
!@ver 1.0 (Q-flux ocean)
IMPLICIT NONE
INTEGER I,J
REAL*8 TGAVE,DWTRO,WTR1O,WTR2O
!@var QTCHNG true if TOCEAN(1) is changed (needed for qflux calculation)
LOGICAL, INTENT(IN) :: QTCHNG
C****
C**** FLAND LAND COVERAGE (1)
C****
C**** TOCEAN 1 OCEAN TEMPERATURE OF FIRST LAYER (C)
C**** 2 MEAN OCEAN TEMPERATURE OF SECOND LAYER (C)
C**** 3 OCEAN TEMPERATURE AT BOTTOM OF SECOND LAYER (C)
C**** RSI RATIO OF OCEAN ICE COVERAGE TO WATER COVERAGE (1)
C****
C**** FWSIM Fresh water sea ice amount (KG/M**2)
C****
C**** RESTRUCTURE OCEAN LAYERS
C****
DO J=1,JM
DO I=1,IMAXJ(J)
IF (FLAND(I,J).GE.1.) CYCLE
IF (Z1OOLD(I,J).GE.Z12O(I,J)) GO TO 140
IF (Z1O(I,J).EQ.Z1OOLD(I,J)) CYCLE
WTR1O=RHOWS*Z1O(I,J)-FWSIM(I,J)
DWTRO=RHOWS*(Z1O(I,J)-Z1OOLD(I,J))
WTR2O=RHOWS*(Z12O(I,J)-Z1O(I,J))
IF (DWTRO.GT.0.) GO TO 120
C**** MIX LAYER DEPTH IS GETTING SHALLOWER
TOCEAN(2,I,J)=TOCEAN(2,I,J)
* +((TOCEAN(2,I,J)-TOCEAN(1,I,J))*DWTRO/WTR2O)
CYCLE
C**** MIX LAYER DEPTH IS GETTING DEEPER
120 IF (QTCHNG) THEN ! change TOCEAN(1)
TGAVE=(TOCEAN(2,I,J)*DWTRO+(2.*TOCEAN(2,I,J)-TOCEAN(3,I,J))
* *WTR2O)/(WTR2O+DWTRO)
TOCEAN(1,I,J)=TOCEAN(1,I,J)+((TGAVE-TOCEAN(1,I,J))
* *DWTRO/WTR1O)
END IF
IF (Z1O(I,J).GE.Z12O(I,J)) GO TO 140
TOCEAN(2,I,J)=TOCEAN(2,I,J)
* +((TOCEAN(3,I,J)-TOCEAN(2,I,J))*DWTRO/(WTR2O+DWTRO))
CYCLE
C**** MIXED LAYER DEPTH IS AT ITS MAXIMUM OR TEMP PROFILE IS UNIFORM
140 TOCEAN(2,I,J)=TOCEAN(1,I,J)
TOCEAN(3,I,J)=TOCEAN(1,I,J)
END DO
END DO
RETURN
END SUBROUTINE OSTRUC
SUBROUTINE OCLIM(end_of_day) 1,19
!@sum OCLIM calculates daily ocean data from ocean/sea ice climatologies
!@auth Original Development Team
!@ver 1.0 (Q-flux ocean or fixed SST)
IMPLICIT NONE
REAL*8, SAVE :: XZO(IM,JM),XZN(IM,JM)
LOGICAL, INTENT(IN) :: end_of_day
INTEGER n,J,I,LSTMON,K,m,m1,JR
REAL*8 RSICSQ,ZIMIN,ZIMAX,Z1OMIN,RSINEW,TIME,FRAC,MSINEW,OPNOCN
* ,TFO
!@var JDLAST julian day that OCLIM was last called
INTEGER, SAVE :: JDLAST=0
!@var IMON0 current month for SST climatology reading
INTEGER, SAVE :: IMON0 = 0
!@var IMON current month for ocean mixed layer climatology reading
INTEGER, SAVE :: IMON = 0
IF (KOCEAN.EQ.1) GO TO 500
if(.not.(end_of_day.or.itime.eq.itimei.or.off_line)) return
C****
C**** OOBS AOST monthly Average Ocean Surface Temperature
C**** BOST 1st order Moment of Ocean Surface Temperature
C**** COST 2nd order Moment of Ocean Surface Temperature
C**** EOST0/1 Ocean Surface Temperature at beginning/end of month
C****
C**** ARSI monthly Average of Ratio of Sea Ice to Water
C**** BRSI 1st order Moment of Ratio of Sea Ice to Water
C**** or .5*((ERSI0-limit)**2+(ERSI1-limit)**2)/(ARSI-limit)
C**** CRSI 2nd order Moment of Ratio of Sea Ice to Water
C**** ERSI0/1 Ratio of Sea Ice to Water at beginning/end of month
C**** KRSI -1: continuous piecewise linear fit at lower limit
C**** 0: quadratic fit
C**** 1: continuous piecewise linear fit at upper limit
C****
C**** CALCULATE DAILY OCEAN DATA FROM CLIMATOLOGY
C****
C**** TOCEAN(1) OCEAN TEMPERATURE (C)
C**** RSI RATIO OF OCEAN ICE COVERAGE TO WATER COVERAGE (1)
C**** MSI OCEAN ICE AMOUNT OF SECOND LAYER (KG/M**2)
C****
C**** READ IN OBSERVED OCEAN DATA
IF (JMON.EQ.IMON0) GO TO 400
IF (IMON0.EQ.0) THEN
C**** READ IN LAST MONTH'S END-OF-MONTH DATA
if (ocn_cycl.ge.1) then
LSTMON=JMON-1
if (lstmon.eq.0) lstmon = 12
call readt (iu_SICE,IM*JM,ERSI0,IM*JM,ERSI0,LSTMON)
if (ocn_cycl.eq.1) then
call readt (iu_OSST,IM*JM,EOST0,IM*JM,EOST0,LSTMON)
else ! if (ocn_cycl.eq.2) then
LSTMON=JMON-1+(JYEAR-IYEAR1)*JMperY
290 read (iu_OSST) M
if (m.lt.lstmon) go to 290
backspace iu_OSST
CALL MREAD (iu_OSST, m,IM*JM,EOST0,IM*JM,EOST0)
WRITE(6,*) 'Read End-of-month ocean data from ',JMON-1,M
IF(M.NE.LSTMON)
& call stop_model('Read error: ocean data',255)
end if
else ! if (ocn_cycl.eq.0) then
LSTMON=JMON-1+(JYEAR-IYEAR1)*JMperY
300 read (iu_OSST) M
if (m.lt.lstmon) go to 300
backspace iu_OSST
310 read (iu_SICE) m
if (m.lt.lstmon) go to 310
backspace iu_SICE
CALL MREAD (iu_OSST, m,IM*JM,EOST0,IM*JM,EOST0)
CALL MREAD (iu_SICE,m1,IM*JM,ERSI0,IM*JM,ERSI0)
WRITE(6,*) 'Read End-of-month ocean data from ',JMON-1,M,M1
IF(M.NE.M1.OR.M.NE.LSTMON)
& call stop_model('Read error: ocean data',255)
end if
ELSE
C**** COPY END-OF-OLD-MONTH DATA TO START-OF-NEW-MONTH DATA
EOST0=EOST1
ERSI0=ERSI1
END IF
C**** READ IN CURRENT MONTHS DATA: MEAN AND END-OF-MONTH
IMON0=JMON
if (ocn_cycl.ge.1) then
if (jmon.eq.1) then
if (ocn_cycl.eq.1) rewind iu_OSST
rewind iu_SICE
read (iu_SICE) ! skip over DM-record
end if
call readt (iu_SICE,0,ARSI,2*im*jm,ARSI,1) ! reads ARSI,ERSI1
if (ocn_cycl.eq.1) then
call readt (iu_OSST,0,AOST,2*im*jm,AOST,1) ! reads AOST,EOST1
else ! if (ocn_cycl.eq.2) then
CALL MREAD (iu_OSST, M,0,AOST,2*IM*JM,AOST) ! READS AOST,EOST1
WRITE(6,*) 'Read in ocean data for month',JMON,M
IF(JMON.NE.MOD(M-1,12)+1)
& call stop_model('Error: Ocean data',255)
end if
else ! if (ocn_cycl.eq.0) then
CALL MREAD (iu_OSST, M,0,AOST,2*IM*JM,AOST) ! READS AOST,EOST1
CALL MREAD (iu_SICE,M1,0,ARSI,2*IM*JM,ARSI) ! READS ARSI,ERSI1
WRITE(6,*) 'Read in ocean data for month',JMON,M,M1
IF(M.NE.M1.OR.JMON.NE.MOD(M-1,12)+1)
& call stop_model('Error: Ocean data',255)
end if
C**** FIND INTERPOLATION COEFFICIENTS (LINEAR/QUADRATIC FIT)
DO J=1,JM
DO I=1,IMAXJ(J)
BOST(I,J)=EOST1(I,J)-EOST0(I,J)
COST(I,J)=3.*(EOST1(I,J)+EOST0(I,J)) - 6.*AOST(I,J)
BRSI(I,J)=0.
CRSI(I,J)=0. ! extreme cases
KRSI(I,J)=0 ! ice constant
IF(ARSI(I,J).LE.0. .or. ARSI(I,J).GE.1.) CYCLE
BRSI(I,J)=ERSI1(I,J)-ERSI0(I,J) ! quadratic fit
CRSI(I,J)=3.*(ERSI1(I,J)+ERSI0(I,J)) - 6.*ARSI(I,J)
IF(ABS(CRSI(I,J)) .GT. ABS(BRSI(I,J))) THEN ! linear fits
RSICSQ=CRSI(I,J)*(ARSI(I,J)*CRSI(I,J) - .25*BRSI(I,J)**2 -
* CRSI(I,J)**2*BY12)
IF(RSICSQ.lt.0.) THEN
C**** RSI uses piecewise linear fit because quadratic fit at apex < 0
KRSI(I,J) = -1
BRSI(I,J) = .5*(ERSI0(I,J)**2 + ERSI1(I,J)**2) / ARSI(I,J)
ELSEIF(RSICSQ.gt.CRSI(I,J)**2) then
C**** RSI uses piecewise linear fit because quadratic fit at apex > 1
KRSI(I,J) = 1
BRSI(I,J) = .5*((ERSI0(I,J)-1.)**2 + (ERSI1(I,J)-1.)**2) /
/ (ARSI(I,J)-1.)
END IF
END IF
END DO
END DO
C**** Calculate OST, RSI and MSI for current day
400 TIME=(JDATE-.5)/(JDendOFM(JMON)-JDendOFM(JMON-1))-.5 ! -.5<TIME<.5
DO J=1,JM
ZIMIN=Z1I+Z2OIM
ZIMAX=2d0
IF(J.GT.JM/2) ZIMAX=3.5d0
DO I=1,IMAXJ(J)
IF (FOCEAN(I,J).gt.0) THEN
C**** OST always uses quadratic fit
TOCEAN(1,I,J)=AOST(I,J)+BOST(I,J)*TIME
* +COST(I,J)*(TIME**2-BY12)
SELECT CASE (KRSI(I,J))
C**** RSI uses piecewise linear fit because quadratic fit at apex < 0
CASE (-1)
IF(ERSI0(I,J)-BRSI(I,J)*(TIME+.5) .gt. 0.) then
RSINEW = ERSI0(I,J) - BRSI(I,J)*(TIME+.5) ! TIME < T0
ELSEIF(ERSI1(I,J)-BRSI(I,J)*(.5-TIME) .gt. 0.) then
RSINEW = ERSI1(I,J) - BRSI(I,J)*(.5-TIME) ! T1 < TIME
ELSE
RSINEW = 0. ! T0 < TIME < T1
END IF
C**** RSI uses piecewise linear fit because quadratic fit at apex > 1
CASE (1)
IF(ERSI0(I,J)-BRSI(I,J)*(TIME+.5) .lt. 1.) then
RSINEW = ERSI0(I,J) - BRSI(I,J)*(TIME+.5) ! TIME < T0
ELSEIF(ERSI1(I,J)-BRSI(I,J)*(.5-TIME) .lt. 1.) then
RSINEW = ERSI1(I,J) - BRSI(I,J)*(.5-TIME) ! T1 < TIME
ELSE
RSINEW = 1. ! T0 < TIME < T1
END IF
C**** RSI uses quadratic fit
CASE (0)
RSINEW=ARSI(I,J)+BRSI(I,J)*TIME+CRSI(I,J)*(TIME**2-BY12)
END SELECT
C**** Set new mass
MSINEW=RHOI*(ZIMIN-Z1I+(ZIMAX-ZIMIN)*RSINEW*DM(I,J))
C**** Ensure that lead fraction is consistent with kocean=1 case
IF (RSINEW.gt.0) THEN
OPNOCN=MIN(0.1d0,FLEADOC*RHOI/(RSINEW*(ACE1I+MSINEW)))
IF (RSINEW.GT.1.-OPNOCN) THEN
RSINEW = 1.-OPNOCN
MSINEW=RHOI*(ZIMIN-Z1I+(ZIMAX-ZIMIN)*RSINEW*DM(I,J))
END IF
END IF
C**** accumulate diagnostics
IF(MSI(I,J).eq.0) MSI(I,J)=MSINEW ! does this happen?
IF (end_of_day.and..not.off_line) THEN
AIJ(I,J,IJ_SMFX)=AIJ(I,J,IJ_SMFX)+
* (SNOWI(I,J)+ACE1I)*(RSINEW-RSI(I,J))+
* RSINEW*MSINEW-RSI(I,J)*MSI(I,J)
AJ(J,J_IMPLM,ITOICE)=AJ(J,J_IMPLM,ITOICE)-FOCEAN(I,J)
* *(MSINEW-MSI(I,J))*RSI(I,J)
AJ(J,J_IMPLH,ITOICE)=AJ(J,J_IMPLH,ITOICE)-FOCEAN(I,J)
* *SUM(HSI(3:4,I,J))*(MSINEW/MSI(I,J)-1.)*RSI(I,J)
AJ(J,J_IMPLM,ITOCEAN)=AJ(J,J_IMPLM,ITOCEAN)-FOCEAN(I,J)
* *(RSINEW-RSI(I,J))*(MSINEW+ACE1I+SNOWI(I,J))
AJ(J,J_IMPLH,ITOCEAN)=AJ(J,J_IMPLH,ITOCEAN)-FOCEAN(I,J)
* *(RSINEW-RSI(I,J))*(SUM(HSI(1:2,I,J))+
* SUM(HSI(3:4,I,J))*(MSINEW/MSI(I,J)))
END IF
C**** adjust enthalpy and salt so temperature/salinity remain constant
HSI(3:4,I,J)=HSI(3:4,I,J)*(MSINEW/MSI(I,J))
SSI(3:4,I,J)=SSI(3:4,I,J)*(MSINEW/MSI(I,J))
RSI(I,J)=RSINEW
MSI(I,J)=MSINEW
C**** WHEN TGO IS NOT DEFINED, MAKE IT A REASONABLE VALUE
TFO=tfrez(sss(i,j))
IF (TOCEAN(1,I,J).LT.TFO) TOCEAN(1,I,J)=TFO
C**** SET DEFAULTS IF NO OCEAN ICE
IF (RSI(I,J).LE.0.) THEN
HSI(1:2,I,J)=(SHI*TFO-LHM*(1.-SSI0))*XSI(1:2)*ACE1I
HSI(3:4,I,J)=(SHI*TFO-LHM*(1.-SSI0))*XSI(3:4)*AC2OIM
SSI(1:2,I,J)=SSI0*XSI(1:2)*ACE1I
SSI(3:4,I,J)=SSI0*XSI(3:4)*AC2OIM
SNOWI(I,J)=0.
GTEMP(1:2,2,I,J)=TFO
END IF
FWSIM(I,J)=RSI(I,J)*(ACE1I+SNOWI(I,J)+MSI(I,J)-SUM(SSI(1:LMI
* ,I,J)))
END IF
END DO
END DO
C**** REPLICATE VALUES AT POLE (for prescribed data only)
IF (FOCEAN(1,JM).gt.0) THEN
DO I=2,IM
SNOWI(I,JM)=SNOWI(1,JM)
TOCEAN(1,I,JM)=TOCEAN(1,1,JM)
RSI(I,JM)=RSI(1,JM)
MSI(I,JM)=MSI(1,JM)
HSI(:,I,JM)=HSI(:,1,JM)
SSI(:,I,JM)=SSI(:,1,JM)
GTEMP(1:2,2,I,JM)=GTEMP(1:2,2,1,JM)
FWSIM(I,JM)=FWSIM(1,JM)
END DO
END IF
RETURN
C****
C**** CALCULATE DAILY OCEAN MIXED LAYER DEPTHS FROM CLIMATOLOGY
C****
C**** SAVE PREVIOUS DAY'S MIXED LAYER DEPTH
500 Z1OOLD=Z1O
C**** COMPUTE Z1O ONLY AT THE END OF A DAY OR AT ItimeI
IF (.not.(end_of_day.or.itime.eq.itimei)) RETURN
C**** Read in climatological ocean mixed layer depths efficiently
IF (JDLAST.eq.0) THEN ! need to read in first month climatology
IMON=1 ! IMON=January
IF (JDAY.LE.16) THEN ! JDAY in Jan 1-15, first month is Dec
CALL READT (iu_OCNML,0,XZO,IM*JM,XZO,12)
REWIND iu_OCNML
ELSE ! JDAY is in Jan 16 to Dec 16, get first month
520 IMON=IMON+1
IF (JDAY.GT.JDmidOFM(IMON) .and. IMON.LE.12) GO TO 520
CALL READT (iu_OCNML,0,XZO,IM*JM,XZO,IMON-1)
IF (IMON.EQ.13) REWIND iu_OCNML
END IF
ELSE ! Do we need to read in second month?
IF (JDAY.NE.JDLAST+1) THEN ! Check that data is read in daily
IF (JDAY.NE.1 .or. JDLAST.NE.365) THEN
WRITE (6,*) 'Incorrect values in OCLIM: JDAY,JDLAST=',JDAY
* ,JDLAST
call stop_model(
& 'ERROR READING IN SETTING OCEAN CLIMATOLOGY',255)
END IF
IMON=IMON-12 ! New year
GO TO 530
END IF
IF (JDAY.LE.JDmidOFM(IMON)) GO TO 530
IMON=IMON+1 ! read in new month of climatological data
XZO = XZN
IF (IMON.EQ.13) REWIND iu_OCNML
END IF
CALL READT (iu_OCNML,0,XZN,IM*JM,XZN,1)
530 JDLAST=JDAY
C**** Interpolate the mixed layer depth z1o to the current day and
C**** limit it to the annual maxmimal mixed layer depth z12o
FRAC = REAL(JDmidOFM(IMON)-JDAY,KIND=8)/
a (JDmidOFM(IMON)-JDmidOFM(IMON-1))
DO J=1,JM
DO I=1,IM
Z1O(I,J)=min( z12o(i,j) , FRAC*XZO(I,J)+(1.-FRAC)*XZN(I,J) )
IF (RSI(I,J)*FOCEAN(I,J).GT.0.and.off_line) THEN
Z1OMIN=1.+FWSIM(I,J)/(RHOWS*RSI(I,J))
IF (Z1OMIN.GT.Z1O(I,J)) THEN
C**** MIXED LAYER DEPTH IS INCREASED TO OCEAN ICE DEPTH + 1 METER
WRITE(6,602) ITime,I,J,JMON,Z1O(I,J),Z1OMIN,z12o(i,j)
602 FORMAT (' INCREASE OF MIXED LAYER DEPTH ',I10,3I4,3F10.3)
Z1O(I,J)=MIN(Z1OMIN, z12o(i,j))
IF (Z1OMIN.GT.Z12O(I,J)) THEN
C**** ICE DEPTH+1>MAX MIXED LAYER DEPTH :
C**** lose the excess mass to the deep ocean
C**** Calculate freshwater mass to be removed, and then any energy/salt
MSINEW=MSI(I,J)*(1.-RHOWS*(Z1OMIN-Z12O(I,J))/(FWSIM(I,J)
* -RSI(I,J)*(ACE1I+SNOWI(I,J)-SUM(SSI(1:2,I,J)))))
HSI(3:4,I,J) = HSI(3:4,I,J)*(MSINEW/MSI(I,J))
SSI(3:4,I,J) = SSI(3:4,I,J)*(MSINEW/MSI(I,J))
AJ(J,J_IMELT,ITOICE)=AJ(J,J_IMPLM,ITOICE)-FOCEAN(I,J)
* *RSI(I,J)*(MSINEW-MSI(I,J))
AJ(J,J_HMELT,ITOICE)=AJ(J,J_IMPLH,ITOICE)-FOCEAN(I,J)
* *RSI(I,J)*SUM(HSI(3:4,I,J))*(MSINEW/MSI(I,J)-1.)
AJ(J,J_SMELT,ITOICE)=AJ(J,J_IMPLH,ITOICE)-FOCEAN(I,J)
* *RSI(I,J)*SUM(SSI(3:4,I,J))*(MSINEW/MSI(I,J)-1.)
AJ(J,J_IMPLM,ITOICE)=AJ(J,J_IMPLM,ITOICE)-FOCEAN(I,J)
* *RSI(I,J)*(MSINEW-MSI(I,J))
AJ(J,J_IMPLH,ITOICE)=AJ(J,J_IMPLH,ITOICE)-FOCEAN(I,J)
* *RSI(I,J)*SUM(HSI(3:4,I,J))*(MSINEW/MSI(I,J)-1.)
JR=JREG(I,J)
AREG(JR,J_IMPLM)=AREG(JR,J_IMPLM)-FOCEAN(I,J)*RSI(I,J)
* *(MSINEW-MSI(I,J))*DXYP(J)
AREG(JR,J_IMPLH)=AREG(JR,J_IMPLH)-FOCEAN(I,J)*RSI(I,J)
* *SUM(HSI(3:4,I,J))*(MSINEW/MSI(I,J)-1.)*DXYP(J)
MSI(I,J)=MSINEW
FWSIM(I,J)=RSI(I,J)*(ACE1I+SNOWI(I,J)+MSI(I,J)-SUM(SSI(1:LMI
* ,I,J)))
END IF
END IF
END IF
END DO
END DO
RETURN
END SUBROUTINE OCLIM
SUBROUTINE OSOURC (ROICE,SMSI,TGW,WTRO,OTDT,RUN0,FODT,FIDT,RVRRUN 1
* ,RVRERUN,EVAPO,EVAPI,TFW,WTRW,RUN4O,ERUN4O
* ,RUN4I,ERUN4I,ENRGFO,ACEFO,ACEFI,ENRGFI)
!@sum OSOURC applies fluxes to ocean in ice-covered and ice-free areas
!@+ ACEFO/I are the freshwater ice amounts,
!@+ ENRGFO/I is total energy (including a salt component)
!@auth Gary Russell
!@ver 1.0
IMPLICIT NONE
!@var TFW freezing temperature for water underlying ice (C)
REAL*8, INTENT(IN) :: TFW
!@var TGW mixed layer temp.(C)
REAL*8, INTENT(INOUT) :: TGW
REAL*8, INTENT(IN) :: ROICE, SMSI, WTRO, EVAPO, EVAPI, RUN0
REAL*8, INTENT(IN) :: FODT, FIDT, OTDT, RVRRUN, RVRERUN
REAL*8, INTENT(OUT) :: ENRGFO, ACEFO, ENRGFI, ACEFI, RUN4O, RUN4I
* , ERUN4O, ERUN4I, WTRW
REAL*8 EIW0, ENRGIW, WTRI1, EFIW, ENRGO0, EOFRZ, ENRGO
REAL*8 WTRW0, ENRGW0, ENRGW, WTRI0, SMSI0
!@var emin min energy deficit required before ice forms (J/m^2)
REAL*8, PARAMETER :: emin=-1d-10
C**** initiallize output
ENRGFO=0. ; ACEFO=0. ; ACEFI=0. ; ENRGFI=0. ; WTRW=WTRO
C**** Calculate previous ice mass (before fluxes applied)
SMSI0=SMSI+RUN0+EVAPI
C**** Calculate extra mass flux to ocean, balanced by deep removal
RUN4O=-EVAPO+RVRRUN ! open ocean
RUN4I=-EVAPI+RVRRUN ! under ice
ERUN4O=RUN4O*TGW*SHW ! corresponding heat flux at bottom (open)
ERUN4I=RUN4I*TGW*SHW ! (under ice)
C**** Calculate heat fluxes to ocean
ENRGO = FODT+OTDT+RVRERUN-ERUN4O ! in open water
ENRGIW = FIDT+OTDT+RVRERUN-ERUN4I ! under ice
C**** Calculate energy in mixed layer (open ocean)
ENRGO0=WTRO*TGW*SHW
EOFRZ=WTRO*TFW*SHW
IF (ENRGO0+ENRGO.GE.EOFRZ+emin) THEN
C**** FLUXES RECOMPUTE TGW WHICH IS ABOVE FREEZING POINT FOR OCEAN
IF (ROICE.LE.0.) THEN
TGW=TGW+ENRGO/(WTRO*SHW)
RETURN
END IF
ELSE
C**** FLUXES COOL TGO TO FREEZING POINT FOR OCEAN AND FORM SOME ICE
ACEFO=(ENRGO0+ENRGO-EOFRZ)/(TFW*(SHI-SHW)-LHM)
ENRGFO=ACEFO*(TFW*SHI/(1.-SSI0)-LHM)
IF (ROICE.LE.0.) THEN
TGW=TFW
RETURN
END IF
END IF
C****
IF (ROICE.le.0) RETURN
C**** Calculate initial energy of mixed layer (before fluxes applied)
WTRI0=WTRO-SMSI0 ! mixed layer mass below ice (kg/m^2)
EIW0=WTRI0*TGW*SHW ! energy of mixed layer below ice (J/m^2)
WTRW0=WTRO-ROICE*SMSI0 ! initial total mixed layer mass
ENRGW0=WTRW0*TGW*SHW ! initial energy in mixed layer
C**** CALCULATE THE ENERGY OF THE WATER BELOW THE ICE AT FREEZING POINT
C**** AND CHECK WHETHER NEW ICE MUST BE FORMED
WTRI1 = WTRO-SMSI ! new mass of ocean (kg/m^2)
EFIW = WTRI1*TFW*SHW ! freezing energy of ocean mass WTRI1
IF (EIW0+ENRGIW .LE. EFIW+emin) THEN ! freezing case
C**** FLUXES WOULD COOL TGW TO FREEZING POINT AND FREEZE SOME MORE ICE
ACEFI = (EIW0+ENRGIW-EFIW)/(TFW*(SHI-SHW)-LHM)
ENRGFI = ACEFI*(TFW*SHI/(1.-SSI0)-LHM) ! energy of frozen ice
END IF
C**** COMBINE OPEN OCEAN AND SEA ICE FRACTIONS TO FORM NEW VARIABLES
WTRW = WTRO -(1.-ROICE)* ACEFO -ROICE*(SMSI+ACEFI)
ENRGW= ENRGW0+(1.-ROICE)*(ENRGO-ENRGFO)+ROICE*(ENRGIW-ENRGFI)
TGW = ENRGW/(WTRW*SHW) ! new ocean temperature
RETURN
END SUBROUTINE OSOURC
END MODULE STATIC_OCEAN
SUBROUTINE init_OCEAN(iniOCEAN,istart) 1,21
!@sum init_OCEAN initiallises ocean variables
!@auth Original Development Team
!@ver 1.0
USE FILEMANAGER
USE PARAM
USE MODEL_COM, only : im,jm,fland,flice,kocean,focean
* ,fearth,iyear1,ioreadnt,jmpery
USE FLUXES, only : gtemp,sss,uosurf,vosurf,uisurf,visurf,ogeoza
USE SEAICE, only : qsfix, osurf_tilt
USE SEAICE_COM, only : snowi
USE STATIC_OCEAN, only : ota,otb,otc,z12o,dm,iu_osst,iu_sice
* ,iu_ocnml,tocean,ocn_cycl,sss0
USE DAGCOM, only : npts,icon_OCE,conpt0
IMPLICIT NONE
LOGICAL :: QCON(NPTS), T=.TRUE. , F=.FALSE.
LOGICAL, INTENT(IN) :: iniOCEAN ! true if starting from ic.
CHARACTER CONPT(NPTS)*10
!@var iu_OHT unit number for reading in ocean heat transports & z12o_max
INTEGER :: iu_OHT,iu_GIC
INTEGER :: I,J,m,ISTART,ioerr
!@var z12o_max maximal mixed layer depth (m) for qflux model
real*8 :: z12o_max
real*8 z1ox(im,jm)
if (kocean.eq.1) then
C**** set conservation diagnostic for ocean heat
CONPT=CONPT0
QCON=(/ F, F, F, T, F, T, F, T, T, F, F/)
CALL SET_CON(QCON,CONPT,"OCN HEAT","(10^6 J/M**2) ",
* "(W/M^2) ",1d-6,1d0,icon_OCE)
end if
if (istart.le.0) return
C**** if starting from AIC/GIC files need additional read for ocean
if (istart.le.2) then
call openunit("GIC",iu_GIC,.true.,.true.)
ioerr=-1
call io_ocean (iu_GIC,ioreadnt,ioerr)
if (ioerr.eq.1) then
WRITE(6,*) "I/O ERROR IN GIC FILE: KUNIT=",iu_GIC
call stop_model("INPUT: GIC READ IN ERROR",255)
end if
call closeunit (iu_GIC)
end if
if (kocean.eq.0) then
call sync_param( "ocn_cycl", ocn_cycl ) ! 1:cycle data 0:dont
C**** set up unit numbers for ocean climatologies
call openunit("OSST",iu_OSST,.true.,.true.)
call openunit("SICE",iu_SICE,.true.,.true.)
if (ocn_cycl.ne.1) then
write(6,*) '********************************************'
write(6,*) '* Make sure that IYEAR1 is consistent with *'
write(6,*) '* the ocean data files OSST and SICE *'
write(6,*) '********************************************'
write(6,*) 'IYEAR1=',IYEAR1
end if
C**** Read in constant factor relating RSI to MSI from sea ice clim.
CALL READT (iu_SICE,0,DM,IM*JM,DM,1)
else ! IF (KOCEAN.eq.1) THEN
C**** DATA FOR QFLUX MIXED LAYER OCEAN RUNS
C**** read in ocean heat transport coefficients
call openunit("OHT",iu_OHT,.true.,.true.)
READ (iu_OHT) OTA,OTB,OTC,z12o_max
WRITE(6,*) "Read ocean heat transports from OHT"
call closeunit (iu_OHT)
C**** Set up unit number of observed mixed layer depth data
call openunit("OCNML",iu_OCNML,.true.,.true.)
C**** find and limit ocean ann max mix layer depths
z12o = 0.
do m=1,jmpery
CALL READT (iu_OCNML,0,z1ox,IM*JM,z1ox,1)
do j=1,jm
do i=1,im
ccc z12o(i,j)=min( z12o_max , max(z12o(i,j),z1ox(i,j)) )
ccc the above line could substitute for next 3 lines w/o any change
if (focean(i,j).gt.0. .and. z1ox(i,j).gt.z12o_max)
* z1ox(i,j)=z12o_max
if (z1ox(i,j).gt.z12o(i,j)) z12o(i,j)=z1ox(i,j)
end do
end do
end do
rewind iu_OCNML
write(6,*) 'Mixed Layer Depths limited to',z12o_max
C**** initialise deep ocean arrays if required
call init_ODEEP(iniOCEAN)
END IF
C**** Set fluxed arrays for oceans
DO J=1,JM
DO I=1,IM
IF (FOCEAN(I,J).gt.0) THEN
GTEMP(1:2,1,I,J)=TOCEAN(1:2,I,J)
SSS(I,J) = SSS0
ELSE
SSS(I,J) = 0.
END IF
C**** For the time being assume zero surface velocities for drag calc
uosurf(i,j)=0. ; vosurf(i,j)=0.
uisurf(i,j)=0. ; visurf(i,j)=0.
C**** Also zero out surface height variations
ogeoza(i,j)=0.
END DO
END DO
C**** keep salinity in sea ice constant for fixed-SST and qflux models
qsfix = .true.
C**** Make sure to use geostrophy for ocean tilt term in ice dynamics
C**** (if required). Since ocean currents are zero, this implies no sea
C**** surface tilt term.
osurf_tilt = 0
C****
RETURN
END SUBROUTINE init_OCEAN
SUBROUTINE daily_OCEAN(end_of_day) 2,8
!@sum daily_OCEAN performs the daily tasks for the ocean module
!@auth Original Development Team
!@ver 1.0
USE CONSTANT, only : twopi,edpery,shw,rhows
USE MODEL_COM, only : im,jm,kocean,focean,jday
USE DAGCOM, only : aij,ij_toc2,ij_tgo2
USE FLUXES, only : gtemp,mlhc,fwsim
USE STATIC_OCEAN, only : tocean,ostruc,oclim,z1o,
* sinang,sn2ang,sn3ang,sn4ang,cosang,cs2ang,cs3ang,cs4ang
IMPLICIT NONE
INTEGER I,J
LOGICAL, INTENT(IN) :: end_of_day
REAL*8 ANGLE
C**** update ocean related climatologies
CALL OCLIM(end_of_day)
C**** Set fourier coefficients for heat transport calculations
IF (KOCEAN.eq.1) THEN
ANGLE=TWOPI*JDAY/EDPERY
SINANG=SIN(ANGLE)
SN2ANG=SIN(2*ANGLE)
SN3ANG=SIN(3*ANGLE)
SN4ANG=SIN(4*ANGLE)
COSANG=COS(ANGLE)
CS2ANG=COS(2*ANGLE)
CS3ANG=COS(3*ANGLE)
CS4ANG=COS(4*ANGLE)
END IF
C**** Only do this at end of the day
IF (KOCEAN.EQ.1.and.end_of_day) THEN
DO J=1,JM
DO I=1,IM
AIJ(I,J,IJ_TOC2)=AIJ(I,J,IJ_TOC2)+TOCEAN(2,I,J)
AIJ(I,J,IJ_TGO2)=AIJ(I,J,IJ_TGO2)+TOCEAN(3,I,J)
END DO
END DO
C**** DO DEEP DIFFUSION IF REQUIRED
CALL ODIFS
C**** RESTRUCTURE THE OCEAN LAYERS
CALL OSTRUC(.TRUE.)
END IF
C**** set gtemp array for ocean temperature
DO J=1,JM
DO I=1,IM
IF (FOCEAN(I,J).gt.0) THEN
GTEMP(1:2,1,I,J) = TOCEAN(1:2,I,J)
MLHC(I,J) = SHW*(Z1O(I,J)*RHOWS-FWSIM(I,J))
END IF
END DO
END DO
C****
RETURN
END SUBROUTINE daily_OCEAN
SUBROUTINE PRECIP_OC 1,8
!@sum PRECIP_OC driver for applying precipitation to ocean fraction
!@auth Original Development Team
!@ver 1.0
USE CONSTANT, only : rhows,shw
USE MODEL_COM, only : im,jm,focean,kocean,itocean,itoice
USE GEOM, only : imaxj,dxyp
USE DAGCOM, only : aj,j_implm,j_implh,oa,areg,jreg
USE FLUXES, only : runpsi,srunpsi,prec,eprec,gtemp,mlhc,melti
* ,emelti,smelti,fwsim
USE SEAICE, only : ace1i
USE SEAICE_COM, only : rsi,msi,snowi
USE STATIC_OCEAN, only : tocean,z1o
IMPLICIT NONE
REAL*8 TGW,PRCP,WTRO,ENRGP,ERUN4,POCEAN,POICE,SNOW
* ,SMSI0,ENRGW,WTRW0,WTRW,RUN0,RUN4,ROICE,SIMELT,ESIMELT
INTEGER I,J,JR
DO J=1,JM
DO I=1,IMAXJ(J)
ROICE=RSI(I,J)
POICE=FOCEAN(I,J)*RSI(I,J)
POCEAN=FOCEAN(I,J)*(1.-RSI(I,J))
JR=JREG(I,J)
IF (FOCEAN(I,J).gt.0) THEN
PRCP=PREC(I,J)
ENRGP=EPREC(I,J)
RUN0=RUNPSI(I,J)-SRUNPSI(I,J) ! fresh water runoff
SIMELT =(MELTI(I,J)-SMELTI(I,J))/(FOCEAN(I,J)*DXYP(J))
ESIMELT=EMELTI(I,J)/(FOCEAN(I,J)*DXYP(J))
OA(I,J,4)=OA(I,J,4)+ENRGP
IF (KOCEAN .EQ. 1) THEN
TGW=TOCEAN(1,I,J)
WTRO=Z1O(I,J)*RHOWS
SNOW=SNOWI(I,J)
SMSI0=FWSIM(I,J)+ROICE*(RUN0-PRCP) ! initial ice
C**** Calculate effect of lateral melt of sea ice
IF (SIMELT.gt.0) THEN
TGW=TGW + (ESIMELT-SIMELT*TGW*SHW)/
* (SHW*(WTRO-SMSI0))
END IF
C**** Additional mass (precip) is balanced by deep removal
RUN4=PRCP
ERUN4=RUN4*TGW*SHW
IF (POICE.LE.0.) THEN
ENRGW=TGW*WTRO*SHW + ENRGP - ERUN4
WTRW =WTRO
ELSE
WTRW0=WTRO -SMSI0
ENRGW=WTRW0*TGW*SHW + (1.-ROICE)*ENRGP - ERUN4
WTRW =WTRW0+ROICE*(RUN0-RUN4)
END IF
TGW=ENRGW/(WTRW*SHW)
TOCEAN(1,I,J)=TGW
AJ(J,J_IMPLM,ITOCEAN)=AJ(J,J_IMPLM,ITOCEAN)+RUN4 *POCEAN
AJ(J,J_IMPLH,ITOCEAN)=AJ(J,J_IMPLH,ITOCEAN)+ERUN4*POCEAN
AJ(J,J_IMPLM,ITOICE) =AJ(J,J_IMPLM,ITOICE) +RUN4 *POICE
AJ(J,J_IMPLH,ITOICE) =AJ(J,J_IMPLH,ITOICE) +ERUN4*POICE
AREG(JR,J_IMPLM)=AREG(JR,J_IMPLM)+RUN4 *FOCEAN(I,J)*DXYP(J)
AREG(JR,J_IMPLH)=AREG(JR,J_IMPLH)+ERUN4*FOCEAN(I,J)*DXYP(J)
MLHC(I,J)=WTRW*SHW ! needed for underice fluxes
END IF
GTEMP(1,1,I,J)=TOCEAN(1,I,J)
END IF
END DO
END DO
RETURN
C****
END SUBROUTINE PRECIP_OC
SUBROUTINE OCEANS 1,10
!@sum OCEANS driver for applying surface fluxes to ocean fraction
!@auth Original Development Team
!@ver 1.0
!@calls OCEAN:OSOURC
USE CONSTANT, only : rhows,shw
USE MODEL_COM, only : im,jm,focean,kocean,jday,dtsrc,itocean
* ,itoice
USE GEOM, only : imaxj,dxyp
USE DAGCOM, only : aj,areg,jreg,j_implm,j_implh,j_oht,oa
USE FLUXES, only : runosi,erunosi,srunosi,e0,e1,evapor,dmsi,dhsi
* ,dssi,flowo,eflowo,gtemp,sss,fwsim,mlhc,gmelt,egmelt
USE SEAICE, only : ace1i,ssi0,tfrez
USE SEAICE_COM, only : rsi,msi,snowi
USE STATIC_OCEAN, only : tocean,z1o,ota,otb,otc,osourc,
* sinang,sn2ang,sn3ang,sn4ang,cosang,cs2ang,cs3ang,cs4ang
IMPLICIT NONE
C**** grid box variables
REAL*8 POCEAN, POICE, DXYPJ, TFO
C**** prognostic variables
REAL*8 TGW, WTRO, SMSI, ROICE
C**** fluxes
REAL*8 EVAPO, EVAPI, FIDT, FODT, OTDT, RVRRUN, RVRERUN, RUN0
C**** output from OSOURC
REAL*8 ERUN4I, ERUN4O, RUN4I, RUN4O, ENRGFO, ACEFO, ACEFI, ENRGFI,
* WTRW
INTEGER I,J,JR
DO J=1,JM
DXYPJ=DXYP(J)
DO I=1,IMAXJ(J)
JR=JREG(I,J)
ROICE=RSI(I,J)
POICE=FOCEAN(I,J)*RSI(I,J)
POCEAN=FOCEAN(I,J)*(1.-RSI(I,J))
IF (FOCEAN(I,J).gt.0) THEN
TGW =TOCEAN(1,I,J)
EVAPO=EVAPOR(I,J,1)
EVAPI=EVAPOR(I,J,2) ! evapor/dew at the ice surface (kg/m^2)
FODT =E0(I,J,1)
SMSI = 0.
IF (ROICE.gt.0) SMSI =FWSIM(I,J)/ROICE
TFO =tfrez(sss(i,j))
C**** get ice-ocean fluxes from sea ice routine
RUN0=RUNOSI(I,J)-SRUNOSI(I,J) ! fw, includes runoff+basal
FIDT=ERUNOSI(I,J)
c SALT=SRUNOSI(I,J)
C**** get river runoff/iceberg melt flux
RVRRUN =( FLOWO(I,J)+ GMELT(I,J))/(FOCEAN(I,J)*DXYPJ)
RVRERUN=(EFLOWO(I,J)+EGMELT(I,J))/(FOCEAN(I,J)*DXYPJ)
OA(I,J,4)=OA(I,J,4)+RVRERUN ! add rvr E to surf. energy budget
IF (KOCEAN .EQ. 1) THEN
WTRO=Z1O(I,J)*RHOWS
OTDT=DTSRC*(OTA(I,J,4)*SN4ANG+OTB(I,J,4)*CS4ANG
* +OTA(I,J,3)*SN3ANG+OTB(I,J,3)*CS3ANG
* +OTA(I,J,2)*SN2ANG+OTB(I,J,2)*CS2ANG
* +OTA(I,J,1)*SINANG+OTB(I,J,1)*COSANG+OTC(I,J))
C**** Calculate the amount of ice formation
CALL OSOURC (ROICE,SMSI,TGW,WTRO,OTDT,RUN0,FODT,FIDT,RVRRUN
* ,RVRERUN,EVAPO,EVAPI,TFO,WTRW,RUN4O,ERUN4O
* ,RUN4I,ERUN4I,ENRGFO,ACEFO,ACEFI,ENRGFI)
C**** Resave prognostic variables
TOCEAN(1,I,J)=TGW
C**** Open Ocean diagnostics
AJ(J,J_OHT ,ITOCEAN)=AJ(J,J_OHT ,ITOCEAN)+OTDT *POCEAN
AJ(J,J_IMPLM,ITOCEAN)=AJ(J,J_IMPLM,ITOCEAN)+RUN4O *POCEAN
AJ(J,J_IMPLH,ITOCEAN)=AJ(J,J_IMPLH,ITOCEAN)+ERUN4O*POCEAN
C**** Ice-covered ocean diagnostics
AJ(J,J_OHT ,ITOICE)=AJ(J,J_OHT ,ITOICE)+OTDT *POICE
AJ(J,J_IMPLM,ITOICE)=AJ(J,J_IMPLM,ITOICE)+RUN4I *POICE
AJ(J,J_IMPLH,ITOICE)=AJ(J,J_IMPLH,ITOICE)+ERUN4I*POICE
C**** regional diagnostics
AREG(JR,J_IMPLM)=AREG(JR,J_IMPLM)+
* (RUN4O *POCEAN+RUN4I *POICE)*DXYPJ
AREG(JR,J_IMPLH)=AREG(JR,J_IMPLH)+
* (ERUN4O*POCEAN+ERUN4I*POICE)*DXYPJ
MLHC(I,J)=SHW*WTRW
ELSE
ACEFO=0 ; ACEFI=0. ; ENRGFO=0. ; ENRGFI=0.
END IF
C**** Store mass and energy fluxes for formation of sea ice
C**** Note ACEFO/I is the freshwater amount of ice.
C**** Add salt for total mass
DMSI(1,I,J)=ACEFO/(1.-SSI0)
DMSI(2,I,J)=ACEFI/(1.-SSI0)
DHSI(1,I,J)=ENRGFO
DHSI(2,I,J)=ENRGFI
DSSI(1,I,J)=SSI0*DMSI(1,I,J)
DSSI(2,I,J)=SSI0*DMSI(2,I,J)
C**** store surface temperatures
GTEMP(1:2,1,I,J)=TOCEAN(1:2,I,J)
END IF
END DO
END DO
RETURN
C****
END SUBROUTINE OCEANS
SUBROUTINE ADVSI_DIAG 1,8
!@sum ADVSI_DIAG adjust diagnostics + mlhc for qflux
!@auth Gavin Schmidt
USE CONSTANT, only : shw,rhows
USE MODEL_COM, only : focean,im,jm,itocean,itoice,kocean,itime
* ,jmon
USE GEOM, only : dxyp,imaxj
USE STATIC_OCEAN, only : tocean,z1o,z12o
USE SEAICE, only : ace1i,lmi
USE SEAICE_COM, only : rsi,msi,hsi,ssi,snowi
USE FLUXES, only : fwsim,msicnv,mlhc
USE DAGCOM, only : aj,areg,J_IMPLM,J_IMPLH,jreg,aij,j_imelt
* ,j_hmelt,j_smelt
IMPLICIT NONE
INTEGER I,J,JR
REAL*8 DXYPJ,RUN4,ERUN4,TGW,POICE,POCEAN,Z1OMIN,MSINEW
IF (KOCEAN.eq.1) THEN ! qflux model
DO J=1,JM
DXYPJ=DXYP(J)
DO I=1,IMAXJ(J)
JR=JREG(I,J)
POICE=FOCEAN(I,J)*RSI(I,J)
POCEAN=FOCEAN(I,J)*(1.-RSI(I,J))
IF (FOCEAN(I,J).gt.0) THEN
TGW = TOCEAN(1,I,J)
RUN4 = MSICNV(I,J)
ERUN4 = TGW*SHW*RUN4
C**** Ensure that we don't run out of ocean if ice gets too thick
IF (POICE.GT.0) THEN
Z1OMIN=1.+FWSIM(I,J)/(RHOWS*RSI(I,J))
IF (Z1OMIN.GT.Z1O(I,J)) THEN
C**** MIXED LAYER DEPTH IS INCREASED TO OCEAN ICE DEPTH + 1 METER
WRITE(6,602) ITime,I,J,JMON,Z1O(I,J),Z1OMIN,z12o(i,j)
602 FORMAT (' INCREASE OF MIXED LAYER DEPTH ',I10,3I4,3F10.3)
Z1O(I,J)=MIN(Z1OMIN, z12o(i,j))
IF (Z1OMIN.GT.Z12O(I,J)) THEN
C**** ICE DEPTH+1>MAX MIXED LAYER DEPTH :
C**** lose the excess mass to the deep ocean
C**** Calculate freshwater mass to be removed, and then any energy/salt
MSINEW=MSI(I,J)*(1.-RHOWS*(Z1OMIN-Z12O(I,J))/(FWSIM(I,J)
* -RSI(I,J)*(ACE1I+SNOWI(I,J)-SUM(SSI(1:2,I,J)))))
HSI(3:4,I,J) = HSI(3:4,I,J)*(MSINEW/MSI(I,J))
SSI(3:4,I,J) = SSI(3:4,I,J)*(MSINEW/MSI(I,J))
AJ(J,J_IMELT,ITOICE)=AJ(J,J_IMPLM,ITOICE)-FOCEAN(I,J)
* *RSI(I,J)*(MSINEW-MSI(I,J))
AJ(J,J_HMELT,ITOICE)=AJ(J,J_IMPLH,ITOICE)-FOCEAN(I,J)
* *RSI(I,J)*SUM(HSI(3:4,I,J))*(MSINEW/MSI(I,J)-1.)
AJ(J,J_SMELT,ITOICE)=AJ(J,J_IMPLH,ITOICE)-FOCEAN(I,J)
* *RSI(I,J)*SUM(SSI(3:4,I,J))*(MSINEW/MSI(I,J)-1.)
AJ(J,J_IMPLM,ITOICE)=AJ(J,J_IMPLM,ITOICE)-FOCEAN(I,J)
* *RSI(I,J)*(MSINEW-MSI(I,J))
AJ(J,J_IMPLH,ITOICE)=AJ(J,J_IMPLH,ITOICE)-FOCEAN(I,J)
* *RSI(I,J)*SUM(HSI(3:4,I,J))*(MSINEW/MSI(I,J)-1.)
AREG(JR,J_IMPLM)=AREG(JR,J_IMPLM)-FOCEAN(I,J)*RSI(I,J)
* *(MSINEW-MSI(I,J))*DXYPJ
AREG(JR,J_IMPLH)=AREG(JR,J_IMPLH)-FOCEAN(I,J)*RSI(I,J)
* *SUM(HSI(3:4,I,J))*(MSINEW/MSI(I,J)-1.)*DXYPJ
MSI(I,J)=MSINEW
FWSIM(I,J)=RSI(I,J)*(ACE1I+SNOWI(I,J)+MSI(I,J)
* -SUM(SSI(1:LMI,I,J)))
END IF
END IF
END IF
MLHC(I,J) = SHW*(Z1O(I,J)*RHOWS-FWSIM(I,J))
C**** Open Ocean diagnostics
AJ(J,J_IMPLM,ITOCEAN)=AJ(J,J_IMPLM,ITOCEAN)+RUN4 *POCEAN
AJ(J,J_IMPLH,ITOCEAN)=AJ(J,J_IMPLH,ITOCEAN)+ERUN4*POCEAN
C**** Ice-covered ocean diagnostics
AJ(J,J_IMPLM,ITOICE)=AJ(J,J_IMPLM,ITOICE)+RUN4 *POICE
AJ(J,J_IMPLH,ITOICE)=AJ(J,J_IMPLH,ITOICE)+ERUN4*POICE
C**** regional diagnostics
AREG(JR,J_IMPLM)=AREG(JR,J_IMPLM)+RUN4 *FOCEAN(I,J)*DXYPJ
AREG(JR,J_IMPLH)=AREG(JR,J_IMPLH)+ERUN4*FOCEAN(I,J)*DXYPJ
END IF
END DO
END DO
END IF
RETURN
END
SUBROUTINE DIAGCO (M) 1,3
!@sum DIAGCO Keeps track of the ocean conservation properties
!@auth Gary Russell/Gavin Schmidt
!@ver 1.0
USE MODEL_COM, only : kocean
USE DAGCOM, only : icon_OCE
IMPLICIT NONE
!@var M index denoting from where DIAGCO is called
INTEGER, INTENT(IN) :: M
C****
C**** THE PARAMETER M INDICATES WHEN DIAGCO IS BEING CALLED
C**** (see DIAGCA)
REAL*8, EXTERNAL :: conserv_OCE
C**** OCEAN POTENTIAL ENTHALPY
IF (KOCEAN.gt.0) CALL conserv_DIAG(M,conserv_OCE,icon_OCE)
C****
RETURN
END SUBROUTINE DIAGCO
SUBROUTINE io_oda(kunit,it,iaction,ioerr) 1,4
!@sum io_oda reads/writes ocean/ice data for initializing deep ocean
!@auth Gavin Schmidt
!@ver 1.0
USE MODEL_COM, only : ioread,iowrite,Itime,im,jm
USE DAGCOM, only : ij_tgo2,aij
USE SEAICE_COM, only : rsi,msi,hsi,ssi
USE STATIC_OCEAN, only : tocean
IMPLICIT NONE
INTEGER kunit !@var kunit unit number of read/write
INTEGER iaction !@var iaction flag for reading or writing to file
!@var IOERR 1 (or -1) if there is (or is not) an error in i/o
INTEGER, INTENT(INOUT) :: IOERR
!@var it input/ouput value of hour
INTEGER, INTENT(INOUT) :: it
INTEGER I,J
SELECT CASE (IACTION)
CASE (:IOWRITE) ! output
WRITE (kunit,err=10) it,((AIJ(I,J,IJ_TGO2),I=1,IM),J=1,JM)
CASE (IOREAD:) ! input
READ (kunit,err=10) it,((AIJ(I,J,IJ_TGO2),I=1,IM),J=1,JM)
END SELECT
RETURN
10 IOERR=1
RETURN
C****
END SUBROUTINE io_oda