module mod_hycom_wavesonly #if defined(USE_ESMF) use ESMF_Mod ! ESMF Framework #endif use mod_xc ! HYCOM communication interface use mod_za ! HYCOM I/O interface use mod_pipe ! HYCOM debugging interface use mod_incupd ! HYCOM incremental update (for data assimilation) use mod_floats ! HYCOM synthetic floats, drifters and moorings use mod_tides ! HYCOM tides use mod_mean ! HYCOM mean archives #if defined(USE_CCSM3) use ccsm3, only : ccsm3_setup_coupling_env ! ccsm3 ocean use ccsm3_exit ! ccsm3 ocean exit use ccsm3_io ! ccsm3 ocean io use ccsm3_forcing ! ccsm3 ocean/cpl6 comms use ccsm3_time_management ! ccsm3 ocean time use ccsm3_global_reductions ! ccsm3 ocean mpi comms use shr_timer_mod ! ccsm3 shared timing routine #endif #if defined(NERSC_VERSION) use mod_year_info, only: year_info,rt,refyear,year_day use mod_forcing_nersc use mod_hycom_nersc use mod_diagnostics use mod_gridp use mod_common_ice use mod_tides_nersc use mod_restart #if defined(WEEKLY_AVERAGE) use mod_average #endif #endif/* NERSC_VERSION*/ c c --- ----------------------------------------- c --- MICOM-based HYbrid Coordinate Ocean Model c --- H Y C O M c --- v e r s i o n 2.2 c --- ----------------------------------------- c implicit none c #if ! defined(USE_CCSM3) include 'common_blocks.h' #endif c #if defined(USE_ESMF) public HYCOM_SetServices #else public HYCOM_Init, HYCOM_Run, HYCOM_Final #endif c logical, save, public :: put_export !set in main program logical, save, public :: get_import !set in main program logical, save, public :: end_of_run !set in HYCOM_Run integer, save, public :: nts_day !set in HYCOM_Init, timesteps/day integer, save, public :: nts_ice !set in HYCOM_Init, timesteps/ice c integer, save, private :: & m,n #if defined(USE_CCSM3) logical, save, private :: & end_month integer, save, private :: & k1m, !k1n in ccsm3_forcing & mm,nn, & termination, !if 1, stop now! & ix0,ix1,ix2 character, save, private :: & ccsm3_string*80 #endif real*8, save, private :: & d1,d2,d3,d4,d2a,d3a,d4a, & ddsurf,ddiagf,dtilef,drstrf,dmeanf, & dske,dskea,dsmr,dsmra,dsms,dsmsa,dsmt,dsmta,dsum,dsuma, & dsumtr(mxtrcr), & dtime,dtime0,dbimon,dmonth,dyear,dyear0, & dsmall,dsmall2 real, save, private :: & smr,sms,smt,sum,smin,smax, tsur, & coord,day1,day2,x,x1,time0,timav,cold,utotp,vtotp real, save, private, allocatable :: & sminy(:),smaxy(:) integer, save, private :: & nstep0,nod, & lt,ma0,ma1,ma2,ma3,mc0,mc1,mc2,mc3, & mk0,mk1,mk2,mk3,mr0,mr1,mr2,mr3,mnth,iofl, #if defined(USE_CCSM3) & jday #else & jday,ihour,iyear #endif logical, save, private :: & linit,diagsv,hisurf,histry,hitile,histmn, & restrt,diag_tide character, save, private :: & intvl*3,c_ydh*14 c real*8 hours1,days1,days6 private hours1,days1,days6 parameter (hours1=1.d0/24.d0,days1=1.d0,days6=6.d0) c c --- tfrz_n = nominal ice melting point (degC) for ice mask real tfrz_n private tfrz_n parameter (tfrz_n=-1.79) !slightly above -1.8 c #if defined(USE_ESMF) c c --- Data types for Import/Export array pointers type ArrayPtrReal2D real(ESMF_KIND_R4), dimension(:,:), pointer :: p end type ArrayPtrReal2D c c --- Attribute names for fields character(ESMF_MAXSTR), save :: & attNameLongName = "long_name", & attNameStdName = "standard_name", & attNameUnits = "units", & attNameSclFac = "scale_factor", & attNameAddOff = "add_offset" c c --- Import Fields integer, parameter :: numImpFields=11 character(ESMF_MAXSTR), save :: impFieldName( numImpFields), & impFieldLongName(numImpFields), & impFieldStdName( numImpFields), & impFieldUnits( numImpFields) real(ESMF_KIND_R4), save :: impFieldSclFac( numImpFields), & impFieldAddOff( numImpFields) c c --- Export Fields integer, parameter :: numExpFields=7 character(ESMF_MAXSTR), save :: expFieldName( numExpFields), & expFieldLongName(numExpFields), & expFieldStdName( numExpFields), & expFieldUnits( numExpFields) real(ESMF_KIND_R4), save :: expFieldSclFac( numExpFields), & expFieldAddOff( numExpFields) c c --- ESMF related variables type(ESMF_Bundle), save :: expBundle, & impBundle type(ESMF_Field), save :: expField(numExpFields), & impField(numImpFields) type(ArrayPtrReal2D), save :: expData( numExpFields), & impData( numImpFields) c type(ESMF_Clock), save :: intClock type(ESMF_VM), save :: vm type(ESMF_DELayout), save :: deLayout integer, save :: petCount, localPet, mpiCommunicator type(ESMF_Grid), save :: grid2D type(ESMF_ArraySpec), save :: arraySpec2Dr real, save, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & sic_import !Sea Ice Concentration &, sitx_import !Sea Ice X-Stress &, sity_import !Sea Ice Y-Stress &, siqs_import !Solar Heat Flux thru Ice to Ocean &, sifh_import !Ice Freezing/Melting Heat Flux &, sifs_import !Ice Freezing/Melting Salt Flux &, sifw_import !Ice Net Water Flux &, sit_import !Sea Ice Temperature &, sih_import !Sea Ice Thickness &, siu_import !Sea Ice X-Velocity &, siv_import !Sea Ice Y-Velocity &, ocn_mask !Ocean Currents Mask logical, save :: & ocn_mask_init #endif contains !!!!#if defined(USE_ESMF) !!!! subroutine HYCOM_SetServices(gridComp, rc) !!!!c !!!! type(ESMF_GridComp) :: gridComp !!!! integer :: rc !!!!c !!!! call ESMF_GridCompSetEntryPoint( !!!! & gridComp, !!!! & ESMF_SETINIT, !!!! & HYCOM_Init, !!!! & ESMF_SINGLEPHASE, !!!! & rc) !!!! call ESMF_GridCompSetEntryPoint( !!!! & gridComp, !!!! & ESMF_SETRUN, !!!! & HYCOM_Run, !!!! & ESMF_SINGLEPHASE, !!!! & rc) !!!! call ESMF_GridCompSetEntryPoint( !!!! & gridComp, !!!! & ESMF_SETFINAL, !!!! & HYCOM_Final, !!!! & ESMF_SINGLEPHASE, !!!! & rc) !!!!c !!!! end subroutine HYCOM_SetServices !!!! !!!! subroutine Setup_ESMF(gridComp, impState, expState, extClock, rc) !!!!c !!!!c --- Calling parameters !!!! type(ESMF_GridComp) :: gridComp !!!! type(ESMF_State) :: impState !!!! type(ESMF_State) :: expState !!!! type(ESMF_Clock) :: extClock !!!! integer :: rc !!!!c !!!!c --- set up ESMF data structures for HYCOM. !!!!c !!!! integer :: i,j !!!! real(ESMF_KIND_R8) :: coord1(itdm),coord2(jtdm) !cell centers !!!! real :: plonrl(itdm),platrl(jtdm) !!!! integer :: cnts(2) !!!! real(ESMF_KIND_R8) :: mgcpd(2), dpd(2) !!!! character(10) :: dimNames(2),dimUnits(2) !!!! type(ESMF_Logical) :: periodic(2) !!!! integer(ESMF_KIND_I4) :: year,month,day,hour,minute !!!! integer(ESMF_KIND_I4) :: sec,msec,usec,nsec !!!! real(8) :: dsec,dmsec,dusec,dnsec !!!! type(ESMF_TimeInterval) :: timeStep, runDuration !!!! type(ESMF_Time) :: startTime !!!! character(ESMF_MAXSTR) :: msg, gridName !!!!c !!!!c --- Report !!!! call ESMF_LogWrite("HYCOM Setup routine called", !!!! & ESMF_LOG_INFO) !!!!!-----call ESMF_LogFlush !!!!c !!!!c Attributes for import fields, identical to CICE export fields !!!! impFieldAddOff(:) = 0.0 !default is no offset !!!! impFieldSclFac(:) = 1.0 !default is no scale factor !!!! !!!! impFieldName( 1) = "sic" !!!! impFieldLongName( 1) = "Sea Ice Concentration" !!!! impFieldStdName( 1) = "sea_ice_area_fraction" !!!! impFieldUnits( 1) = "1" !!!! impFieldName( 2) = "sitx" !!!! impFieldLongName( 2) = "Sea Ice X-Stress" !!!! impFieldStdName( 2) = "downward_x_stress_at_sea_ice_base" !!!! impFieldSclFac( 2) = -1.0 !field is upward !!!! impFieldUnits( 2) = "Pa" !!!! impFieldName( 3) = "sity" !!!! impFieldLongName( 3) = "Sea Ice Y-Stress" !!!! impFieldStdName( 3) = "downward_y_stress_at_sea_ice_base" !!!! impFieldSclFac( 3) = -1.0 !field is upward !!!! impFieldUnits( 3) = "Pa" !!!! impFieldName( 4) = "siqs" !!!! impFieldLongName( 4) = "Solar Heat Flux thru Ice to Ocean" !!!! impFieldStdName( 4) = "downward_sea_ice_basal_solar_heat_flux" !!!! impFieldUnits( 4) = "W m-2" !!!! impFieldName( 5) = "sifh" !!!! impFieldLongName( 5) = "Ice Freezing/Melting Heat Flux" !!!! impFieldStdName( 5) = "upward_sea_ice_basal_heat_flux" !!!! impFieldSclFac( 5) = -1.0 !field is downward !!!! impFieldUnits( 5) = "W m-2" !!!! impFieldName( 6) = "sifs" !!!! impFieldLongName( 6) = "Ice Freezing/Melting Salt Flux" !!!! impFieldStdName( 6) = "downward_sea_ice_basal_salt_flux" !!!! impFieldUnits( 6) = "kg m-2 s-1" !!!! impFieldName( 7) = "sifw" !!!! impFieldLongName( 7) = "Ice Net Water Flux" !!!! impFieldStdName( 7) = "downward_sea_ice_basal_water_flux" !!!! impFieldUnits( 7) = "kg m-2 s-1" !!!! impFieldName( 8) = "sit" !diagnostic !!!! impFieldLongName( 8) = "Sea Ice Temperature" !!!! impFieldStdName( 8) = "sea_ice_temperature" !!!! impFieldAddOff( 8) = +273.15 !field is in degC !!!! impFieldUnits( 8) = "K" !!!! impFieldName( 9) = "sih" !diagnostic !!!! impFieldLongName( 9) = "Sea Ice Thickness" !!!! impFieldStdName( 9) = "sea_ice_thickness" !!!! impFieldUnits( 9) = "m" !!!! impFieldName( 10) = "siu" !diagnostic !!!! impFieldLongName(10) = "Sea Ice X-Velocity" !!!! impFieldStdName( 10) = "sea_ice_x_velocity" !!!! impFieldUnits( 10) = "m s-1" !!!! impFieldName( 11) = "siv" !diagnostic !!!! impFieldLongName(11) = "Sea Ice Y-Velocity" !!!! impFieldStdName( 11) = "sea_ice_y_velocity" !!!! impFieldUnits( 11) = "m s-1" !!!!* !!!!* impFieldName( 12) = "patm" !!!!* impFieldLongName(12) = "Surface Air Pressure" !!!!* impFieldStdName( 12) = "surface_air_pressure" !!!!* impFieldUnits( 12) = "Pa" !!!!* impFieldName( 13) = "xwnd" !!!!* impFieldLongName(13) = "X-Wind" !!!!* impFieldStdName( 13) = "x_wind" !!!!* impFieldUnits( 13) = "m s-1" !!!!* impFieldName( 14) = "ywnd" !!!!* impFieldLongName(14) = "Y-Wind" !!!!* impFieldStdName( 14) = "y_wind" !!!!* impFieldUnits( 14) = "m s-1" !!!!c !!!!c Attributes for export fields, identical to CICE import fields !!!! expFieldAddOff(:) = 0.0 !default is no offset !!!! expFieldSclFac(:) = 1.0 !default is no scale factor !!!! !!!! expFieldName( 1) = "sst" !!!! expFieldLongName( 1) = "Sea Surface Temperature" !!!! expFieldStdName( 1) = "sea_surface_temperature" !!!! expFieldAddOff( 1) = +273.15 !field is in degC !!!! expFieldUnits( 1) = "K" !!!! expFieldName( 2) = "sss" !!!! expFieldLongName( 2) = "Sea Surface Salinity" !!!! expFieldStdName( 2) = "sea_surface_salinity" !!!! expFieldUnits( 2) = "1e-3" !!!! expFieldName( 3) = "ssu" !!!! expFieldLongName( 3) = "Sea Surface X-Current" !!!! expFieldStdName( 3) = "sea_water_x_velocity" !!!! expFieldUnits( 3) = "m s-1" !!!! expFieldName( 4) = "ssv" !!!! expFieldLongName( 4) = "Sea Surface Y-Current" !!!! expFieldStdName( 4) = "sea_water_y_velocity" !!!! expFieldUnits( 4) = "m s-1" !!!! expFieldName( 5) = "ssh" !!!! expFieldLongName( 5) = "Sea Surface Height" !!!! expFieldStdName( 5) = "sea_surface_height_above_sea_level" !!!! expFieldUnits( 5) = "m" !!!! expFieldName( 6) = "ssfi" !!!! expFieldLongName( 6) = "Oceanic Heat Flux Available to Sea Ice" !!!! expFieldStdName( 6) = "upward_sea_ice_basal_available_heat_flux" !!!! expFieldSclFac( 6) = -1.0 !field is downward !!!! expFieldUnits( 6) = "W m-2" !!!! expFieldName( 7) = "mlt" !diagnostic !!!! expFieldLongName( 7) = "Ocean Mixed Layer Thickness" !!!! expFieldStdName( 7) = "ocean_mixed_layer_thickness" !!!! expFieldUnits( 7) = "m" !!!!c !!!!c Create a DE layout to match HYCOM layout !!!! deLayout = ESMF_DELayoutCreate(vm, !!!! & deCountList=(/ipr, jpr/), !!!! & rc=rc) !!!! if (ESMF_LogMsgFoundError(rc, !!!! & "Setup_ESMF: DELayoutCreate failed", rc)) !!!! & call ESMF_Finalize(rc=rc) !!!!c !!!!c Create array specifications !!!! call ESMF_ArraySpecSet(arraySpec2Dr, !!!! & rank=2, !!!! & type=ESMF_DATA_REAL, !!!! & kind=ESMF_R4, !!!! & rc=rc) !!!! if (ESMF_LogMsgFoundError(rc, !!!! & "Setup_ESMF: ArraySpecSet failed", rc)) !!!! & call ESMF_Finalize(rc=rc) !!!!c !!!!c Create an ESMF grid that matches the HYCOM 2D grid !!!!#if defined(ESMF_CURVILINEAR) !!!!c Use indices as cartesian coordinates !!!!c !!!! cnts( 1)=itdm; cnts( 2)=jtdm; !!!! mgcpd( 1)=1.d0; mgcpd( 2)=1.d0; !!!! dpd( 1)=1.d0; dpd( 2)=1.d0; !!!! dimNames(1)="index_i"; dimNames(2)="index_j"; !!!! dimUnits(1)="1"; dimUnits(2)="1"; !!!! periodic(1)=ESMF_FALSE; periodic(2)=ESMF_FALSE; !!!! gridName="HYCOM X-Y Cell-Center Grid" !!!! grid2D = ESMF_GridCreateHorzXYUni( !!!! & counts =cnts, !!!! & minGlobalCoordPerDim=mgcpd, !!!! & deltaPerDim=dpd, !!!! & horzStagger=ESMF_GRID_HORZ_STAGGER_A, !CELL_CENTER only !!!! & dimNames =dimNames, !!!! & dimUnits =dimUnits, !!!! & coordOrder =ESMF_COORD_ORDER_XYZ, !!!! & periodic =periodic, !!!! & name =trim(gridName), !!!! & rc=rc) !!!! if (ESMF_LogMsgFoundError(rc, !!!! & "Setup_ESMF: GridCreateHorzXYUni", rc)) !!!! & call ESMF_Finalize(rc=rc) !!!!#else !!!!c Rectilinear, lat-lon, coordinates !!!! call xclget(plonrl,itdm, plon,1,1,1,0, 0) !!!! call xclget(platrl,jtdm, plat,1,1,0,1, 0) !!!! coord1(1:itdm) = plonrl(1:itdm) !!!! coord2(1:jtdm) = platrl(1:jtdm) !!!! dimNames(1)="longitude"; dimNames(2)="latitude"; !!!! dimUnits(1)="degrees_east"; dimUnits(2)="degrees_north"; !!!! if (plonrl(itdm)-plonrl(1) .gt. 350.0) then !probably global !!!! periodic(1)=ESMF_TRUE !!!! periodic(2)=ESMF_FALSE !!!! else !!!! periodic(1)=ESMF_FALSE !!!! periodic(2)=ESMF_FALSE !!!! endif !!!! gridName="HYCOM Lat-Lon Cell-Center Grid" !!!! grid2D = ESMF_GridCreateHorzLatLon( !!!! & coord1 =coord1(1:itdm), !!!! & coord2 =coord2(1:jtdm), !!!! & horzStagger=ESMF_GRID_HORZ_STAGGER_A, !CELL_CENTER only !!!! & dimNames =dimNames, !!!! & dimUnits =dimUnits, !!!! & coordOrder =ESMF_COORD_ORDER_XYZ, !!!! & periodic =periodic, !!!! & name =trim(gridName), !!!! & rc=rc) !!!! if (ESMF_LogMsgFoundError(rc, !!!! & "Setup_ESMF: GridCreateHorzLatLon", rc)) !!!! & call ESMF_Finalize(rc=rc) !!!!#endif !!!!c !!!!c Distribute the grid on the DE layout using HYCOM distribution !!!!c Only Separable 2-D block distributions with no land skipping are allowed !!!! call ESMF_GridDistribute(grid2D, !!!! & deLayout =deLayout, !!!! & countsPerDEDim1=countde1(1:ipr), !!!! & countsPerDEDim2=countde2(1:jpr), !!!! & rc=rc) !!!! if (ESMF_LogMsgFoundError(rc, !!!! & "Setup_ESMF: GridDistribute", rc)) !!!! & call ESMF_Finalize(rc=rc) !!!!c !!!!c Associate grid with ESMF gridded component !!!! call ESMF_GridCompSet(gridComp, grid=grid2D, rc=rc) !!!! if (ESMF_LogMsgFoundError(rc, !!!! & "Setup_ESMF: GridCompSet", rc)) !!!! & call ESMF_Finalize(rc=rc) !!!!c Setup export fields, bundles & state !!!! do i = 1,numExpFields !!!! expField(i) = ESMF_FieldCreate(grid2D, arraySpec2Dr, !!!! & allocFlag =ESMF_ALLOC, !!!! & horzRelloc=ESMF_CELL_CENTER, !!!! & haloWidth =0, !!!! & name =trim(expFieldName(i)), !!!! & rc=rc) !!!! call ESMF_FieldSetAttribute(expField(i), !!!! & trim(attNameLongName), trim(expFieldLongName(i)), rc=rc) !!!! call ESMF_FieldSetAttribute(expField(i), !!!! & trim(attNameStdName), trim(expFieldStdName(i)), rc=rc) !!!! call ESMF_FieldSetAttribute(expField(i), !!!! & trim(attNameUnits), trim(expFieldUnits(i)), rc=rc) !!!! call ESMF_FieldSetAttribute(expField(i), !!!! & trim(attNameAddOff), expFieldAddOff(i), rc=rc) !!!! call ESMF_FieldSetAttribute(expField(i), !!!! & trim(attNameSclFac), expFieldSclFac(i), rc=rc) !!!! call ESMF_FieldGetDataPointer(expField(i), !!!! & expData(i)%p, copyFlag=ESMF_DATA_REF, rc=rc) !!!! expData(i)%p(:,:) = 0.0 !!!! enddo !!!!c !!!!c Create bundle from list of fields !!!! expBundle = ESMF_BundleCreate(numExpFields, !!!! & expField(:), name="HYCOM", !!!! & rc=rc) !!!!c !!!!c Add bundle to the export state !!!! call ESMF_StateAddBundle(expState, expBundle, rc=rc) !!!!c !!!!c Setup import fields, bundles & state !!!! do i = 1,numImpFields !!!! impField(i) = ESMF_FieldCreate(grid2D, arraySpec2Dr, !!!! & allocFlag =ESMF_ALLOC, !!!! & horzRelloc=ESMF_CELL_CENTER, !!!! & haloWidth =0, !!!! & name =trim(impFieldName(i)), !!!! & rc=rc) !!!! call ESMF_FieldSetAttribute(impField(i), !!!! & trim(attNameLongName), trim(impFieldLongName(i)), rc=rc) !!!! call ESMF_FieldSetAttribute(impField(i), !!!! & trim(attNameStdName), trim(impFieldStdName(i)), rc=rc) !!!! call ESMF_FieldSetAttribute(impField(i), !!!! & trim(attNameUnits), trim(impFieldUnits(i)), rc=rc) !!!! call ESMF_FieldSetAttribute(impField(i), !!!! & trim(attNameAddOff), impFieldAddOff(i), rc=rc) !!!! call ESMF_FieldSetAttribute(impField(i), !!!! & trim(attNameSclFac), impFieldSclFac(i), rc=rc) !!!! call ESMF_FieldGetDataPointer(impField(i), !!!! & impData(i)%p, copyFlag=ESMF_DATA_REF, rc=rc) !!!! impData(i)%p(:,:) = 0.0 !!!! enddo !!!! sic_import(:,:) = 0.0 !Sea Ice Concentration !!!! sitx_import(:,:) = 0.0 !Sea Ice X-Stress !!!! sity_import(:,:) = 0.0 !Sea Ice Y-Stress !!!! siqs_import(:,:) = 0.0 !Solar Heat Flux thru Ice to Ocean !!!! sifh_import(:,:) = 0.0 !Ice Freezing/Melting Heat Flux !!!! sifs_import(:,:) = 0.0 !Ice Freezing/Melting Salt Flux !!!! sifw_import(:,:) = 0.0 !Ice Net Water Flux !!!! sit_import(:,:) = 0.0 !Sea Ice Temperature !!!! sih_import(:,:) = 0.0 !Sea Ice Thickness !!!! siu_import(:,:) = 0.0 !Sea Ice X-Velocity !!!! siv_import(:,:) = 0.0 !Sea Ice Y-Velocity !!!!c !!!!c Create bundle from list of fields !!!! impBundle = ESMF_BundleCreate(numImpFields, !!!! & impField(:), name="HYCOM", !!!! & rc=rc) !!!!c !!!!c Add bundle to the import state !!!! call ESMF_StateAddBundle(impState, impBundle, rc=rc) !!!!c !!!! ocn_mask_init = .true. !still need to initialize ocn_mask !!!!c !!!! end subroutine Setup_ESMF !!!! !!!! subroutine Export_ESMF !!!!c !!!!c --- Fill export state. !!!!c --- Calculate ssfi "in place" !!!!c !!!! integer i,j,k !!!! real ssh2m !!!! real tmxl,smxl,umxl,vmxl,hfrz,tfrz,t2f,ssfi !!!! real dp1,usur1,vsur1,psur1,dp2,usur2,vsur2,psur2 !!!!c !!!!c --- Report !!!! call ESMF_LogWrite("HYCOM Export routine called", !!!! & ESMF_LOG_INFO) !!!!!-----call ESMF_LogFlush !!!!c !!!! margin = 0 !!!!c !!!! if (ocn_mask_init) then !very 1st call to this routine !!!! ocn_mask_init = .false. !!!!c !!!! if (iceflg.eq.4) then !!!! ocn_mask(:,:) = 0.0 !export ocean currents nowhere !!!! elseif (nestfq.ne.0.0) then !!!!c export ocean currents away from open boundaries !!!! do j= 1,jj !!!! do i= 1,ii !!!! if (rmunv(i,j).ne.0.0) then !!!! ocn_mask(i,j) = 0.0 !!!! else !!!! ocn_mask(i,j) = 1.0 !!!! endif !!!! enddo !i !!!! enddo !j !!!! do i= 1,10 !!!! call psmooth(ocn_mask,0) !not efficient, but only done once !!!! enddo !i !!!! else !!!! ocn_mask(:,:) = 1.0 !export ocean currents everywhere !!!! endif !!!! endif !ocn_mask_init !!!!c !!!!c --- Assume Export State is as defined in Setup_ESMF !!!!c --- Average two time levels since (the coupling frequency) icefrq >> 2 !!!!c !!!! ssh2m = 1.0/g !!!! do j= 1,jj !!!! do i= 1,ii !!!! if (ip(i,j).eq.1) then !!!!c --- quantities for available freeze/melt heat flux !!!!c --- relax to tfrz with e-folding time of icefrq time steps !!!!c --- assuming the effective surface layer thickness is hfrz !!!!c --- multiply by dpbl(i,j)/hfrz to get the actual e-folding time !!!! hfrz = min( thkfrz*onem, dpbl(i,j) ) !!!! t2f = (spcifh*hfrz)/(baclin*icefrq*g) !!!! smxl = 0.5*(saln(i,j,1,n)+saln(i,j,1,m)) !!!! tmxl = 0.5*(temp(i,j,1,n)+temp(i,j,1,m)) !!!! tfrz = tfrz_0 + smxl*tfrz_s !salinity dependent freezing point !!!! ssfi = (tfrz-tmxl)*t2f !W/m^2 into ocean !!!!c --- average currents over top 10m !!!! usur1 = 0.0 !!!! vsur1 = 0.0 !!!! psur1 = 0.0 !!!! usur2 = 0.0 !!!! vsur2 = 0.0 !!!! psur2 = 0.0 !!!! do k= 1,kk !!!! dp1 = min( dp(i,j,k,1), max( 0.0, tenm - psur1 ) ) !!!! usur1 = usur1 + dp1*(u(i,j,k,1)+u(i+1,j,k,1)) !!!! vsur1 = vsur1 + dp1*(v(i,j,k,1)+v(i,j+1,k,1)) !!!! psur1 = psur1 + dp1 !!!! dp2 = min( dp(i,j,k,2), max( 0.0, tenm - psur2 ) ) !!!! usur2 = usur2 + dp2*(u(i,j,k,2)+u(i+1,j,k,2)) !!!! vsur2 = vsur2 + dp2*(v(i,j,k,2)+v(i,j+1,k,2)) !!!! psur2 = psur2 + dp2 !!!! if (max(dp1,dp2).eq.0.0) then !!!! exit !!!! endif !!!! enddo !!!! umxl = 0.25*( usur1/psur1 + ubavg(i, j,1) + !!!! + ubavg(i+1,j,1) + !!!! & usur2/psur2 + ubavg(i, j,2) + !!!! & ubavg(i+1,j,2) ) !!!! vmxl = 0.25*( vsur1/psur1 + vbavg(i,j, 1) + !!!! & vbavg(i,j+1,1) + !!!! & vsur2/psur2 + vbavg(i,j, 2) + !!!! & vbavg(i,j+1,2) ) !!!! expData(1)%p(i,j) = tmxl !!!! expData(2)%p(i,j) = smxl !!!! expData(3)%p(i,j) = umxl*ocn_mask(i,j) !!!! expData(4)%p(i,j) = vmxl*ocn_mask(i,j) !!!! expData(5)%p(i,j) = ssh2m*srfhgt(i,j) !ssh in m !!!! expData(6)%p(i,j) = max(-1000.0,min(1000.0,ssfi)) !as in CICE !!!! expData(7)%p(i,j) = dpbl(i,j)*qonem !!!! endif !ip !!!! enddo !i !!!! enddo !j !!!!c !!!! end subroutine Export_ESMF !!!! !!!! subroutine Import_ESMF !!!!c !!!!c --- Extract import state. !!!!c !!!! integer i,j !!!!c !!!!c --- Report !!!! call ESMF_LogWrite("HYCOM Import routine called", !!!! & ESMF_LOG_INFO) !!!!!-----call ESMF_LogFlush !!!!c !!!!c --- Assume Import State is as defined in Setup_ESMF !!!!c !!!! do j= 1,jj !!!! do i= 1,ii !!!! if (ip(i,j).eq.1) then !!!! sic_import(i,j) = impData( 1)%p(i,j) !Sea Ice Concentration !!!! sitx_import(i,j) = impData( 2)%p(i,j) !Sea Ice X-Stress !!!! sity_import(i,j) = impData( 3)%p(i,j) !Sea Ice Y-Stress !!!! siqs_import(i,j) = impData( 4)%p(i,j) !Solar Heat Flux thru Ice to Ocean !!!! sifh_import(i,j) = impData( 5)%p(i,j) !Ice Freezing/Melting Heat Flux !!!! sifs_import(i,j) = impData( 6)%p(i,j) !Ice Freezing/Melting Salt Flux !!!! sifw_import(i,j) = impData( 7)%p(i,j) !Ice Net Water Flux !!!! sit_import(i,j) = impData( 8)%p(i,j) !Sea Ice Temperature !!!! sih_import(i,j) = impData( 9)%p(i,j) !Sea Ice Thickness !!!! siu_import(i,j) = impData(10)%p(i,j) !Sea Ice X-Velocity !!!! siv_import(i,j) = impData(11)%p(i,j) !Sea Ice Y-Velocity !!!! if (iceflg.ge.2 .and. icmflg.ne.3) then !!!! covice(i,j) = impData(1)%p(i,j) !Sea Ice Concentration !!!! si_c(i,j) = impData(1)%p(i,j) !Sea Ice Concentration !!!! if (covice(i,j).gt.0.0) then !!!! si_tx(i,j) = -impData( 2)%p(i,j) !Sea Ice X-Stress into ocean !!!! si_ty(i,j) = -impData( 3)%p(i,j) !Sea Ice Y-Stress into ocean !!!! fswice(i,j) = impData( 4)%p(i,j) !Solar Heat Flux thru Ice to Ocean !!!! flxice(i,j) = fswice(i,j) + !!!! & impData( 5)%p(i,j) !Ice Freezing/Melting Heat Flux !!!! sflice(i,j) = impData( 6)%p(i,j)*1.e3 - !!!! & impData( 7)%p(i,j)*saln(i,j,1,n) !!!! !Ice Virtual Salt Flux !!!! temice(i,j) = impData( 8)%p(i,j) !Sea Ice Temperature !!!! thkice(i,j) = impData( 9)%p(i,j) !Sea Ice Thickness !!!! si_u(i,j) = impData(10)%p(i,j) !Sea Ice X-Velocity !!!! si_v(i,j) = impData(11)%p(i,j) !Sea Ice Y-Velocity !!!! else !!!! si_tx(i,j) = 0.0 !!!! si_ty(i,j) = 0.0 !!!! fswice(i,j) = 0.0 !!!! flxice(i,j) = 0.0 !!!! sflice(i,j) = 0.0 !!!! temice(i,j) = 0.0 !!!! thkice(i,j) = 0.0 !!!! si_u(i,j) = 0.0 !!!! si_v(i,j) = 0.0 !!!! endif !covice !!!! elseif (iceflg.ge.2 .and. icmflg.eq.3) then !!!! si_c(i,j) = impData( 1)%p(i,j) !Sea Ice Concentration !!!! if (si_c(i,j).gt.0.0) then !!!! si_tx(i,j) = -impData( 2)%p(i,j) !Sea Ice X-Stress into ocean !!!! si_ty(i,j) = -impData( 3)%p(i,j) !Sea Ice Y-Stress into ocean !!!! si_h(i,j) = impData( 9)%p(i,j) !Sea Ice Thickness !!!! si_t(i,j) = impData( 8)%p(i,j) !Sea Ice Temperature !!!! si_u(i,j) = impData(10)%p(i,j) !Sea Ice X-Velocity !!!! si_v(i,j) = impData(11)%p(i,j) !Sea Ice Y-Velocity !!!! else !!!! si_tx(i,j) = 0.0 !!!! si_ty(i,j) = 0.0 !!!! si_h(i,j) = 0.0 !!!! si_t(i,j) = 0.0 !!!! si_u(i,j) = 0.0 !!!! si_v(i,j) = 0.0 !!!! endif !covice !!!! endif !iceflg>=2 (icmflg) !!!! endif !ip !!!! enddo !i !!!! enddo !j !!!!c !!!! end subroutine Import_ESMF !!!! !!!! subroutine Archive_ESMF(iyear,jday,ihour) !!!! integer iyear,jday,ihour !!!!c !!!!c --- Create a HYCOM "archive-like" file from Import/Export state.. !!!!c --- Import state may not be at the same time as Export. !!!!c !!!! character*8 cname !!!! character*80 cfile !!!! integer i,j,k,nop,nopa !!!! real coord,xmin,xmax !!!!c !!!! write(cfile,'(a,i4.4,a1,i3.3,a1,i2.2)') !!!! & 'arche.',iyear,'_',jday,'_',ihour !!!! nopa=13 !!!! nop =13+uoff !!!!c !!!! call zaiopf(trim(cfile)//'.a', 'new', nopa) !!!! if (mnproc.eq.1) then !!!! open (unit=nop,file=trim(cfile)//'.b',status='new') !uoff+13 !!!! write(nop,116) ctitle,iversn,iexpt,yrflag,itdm,jtdm !!!! call flush(nop) !!!! endif !1st tile !!!! 116 format (a80/a80/a80/a80/ !!!! & i5,4x,'''iversn'' = hycom version number x10'/ !!!! & i5,4x,'''iexpt '' = experiment number x10'/ !!!! & i5,4x,'''yrflag'' = days in year flag'/ !!!! & i5,4x,'''idm '' = longitudinal array size'/ !!!! & i5,4x,'''jdm '' = latitudinal array size'/ !!!! & 'field time step model day', !!!! & ' k dens min max') !!!!c !!!!c --- surface fields !!!!c !!!! coord=0.0 !!!! do k= 1,numExpFields !!!! do j= 1,jj !!!! do i= 1,ii !!!! if (ip(i,j).eq.1) then !!!! util1(i,j) = expData(k)%p(i,j) !!!! endif !ip !!!! enddo !i !!!! enddo !j !!!! cname = expFieldName(k)(1:8) !!!! call zaiowr(util1,ip,.true., !!!! & xmin,xmax, nopa, .false.) !!!! if (mnproc.eq.1) then !!!! write (nop,117) cname,nstep,time,0,coord,xmin,xmax !!!! call flush(nop) !!!! endif !1st tile !!!! enddo !k !!!! do j= 1,jj !!!! do i= 1,ii !!!! if (ip(i,j).eq.1) then !!!! util2(i,j) = impData(1)%p(i,j) !ice concentration !!!! else !!!! util2(i,j) = 0.0 !!!! endif !ip !!!! enddo !i !!!! enddo !j !!!! cname = impFieldName(1)(1:8) !!!! call zaiowr(util2,ip,.true., !!!! & xmin,xmax, nopa, .false.) !!!! if (mnproc.eq.1) then !!!! write (nop,117) cname,nstep,time,0,coord,xmin,xmax !!!! call flush(nop) !!!! endif !1st tile !!!! do k= 2,3 !si_tx,si_ty !!!! do j= 1,jj !!!! do i= 1,ii !!!! if (util2(i,j).ne.0.0) then !!!! util1(i,j) = -impData(k)%p(i,j) !into ocean !!!! else !!!! util1(i,j) = 0.0 !!!! endif !ice:no-ice !!!! enddo !i !!!! enddo !j !!!! cname = impFieldName(k)(1:8) !!!! call zaiowr(util1,ip,.true., !!!! & xmin,xmax, nopa, .false.) !!!! if (mnproc.eq.1) then !!!! write (nop,117) cname(1:4)//'down', !!!! & nstep,time,0,coord,xmin,xmax !!!! call flush(nop) !!!! endif !1st tile !!!! enddo !k !!!! do k= 4,7 !fluxes !!!! do j= 1,jj !!!! do i= 1,ii !!!! if (util2(i,j).ne.0.0) then !!!! util1(i,j) = util2(i,j)*impData(k)%p(i,j) !!!! else !!!! util1(i,j) = huge !mask where there is no ice !!!! endif !ice:no-ice !!!! enddo !i !!!! enddo !j !!!! cname = impFieldName(k)(1:8) !!!! call zaiowr(util1,ip,.false., !mask on ice !!!! & xmin,xmax, nopa, .false.) !!!! if (mnproc.eq.1) then !!!! write (nop,117) cname,nstep,time,0,coord,xmin,xmax !!!! call flush(nop) !!!! endif !1st tile !!!! enddo !k !!!! do k= 8,numImpFields !!!! do j= 1,jj !!!! do i= 1,ii !!!! if (util2(i,j).ne.0.0) then !!!! util1(i,j) = impData(k)%p(i,j) !!!! else !!!! util1(i,j) = huge !mask where there is no ice !!!! endif !ice:no-ice !!!! enddo !i !!!! enddo !j !!!! cname = impFieldName(k)(1:8) !!!! call zaiowr(util1,ip,.false., !mask on ice !!!! & xmin,xmax, nopa, .false.) !!!! if (mnproc.eq.1) then !!!! write (nop,117) cname,nstep,time,0,coord,xmin,xmax !!!! call flush(nop) !!!! endif !1st tile !!!! enddo !k !!!! do j= 1,jj !!!! do i= 1,ii !!!! if (util2(i,j).ne.0.0) then !!!! util1(i,j) = impData( 6)%p(i,j)*1.e3 - !!!! & impData( 7)%p(i,j)*saln(i,j,1,n) !virtual salt flux !!!! else !!!! util1(i,j) = huge !mask where there is no ice !!!! endif !ice:no-ice !!!! enddo !i !!!! enddo !j !!!! cname = 'surtx ' !!!! call zaiowr(surtx,ip,.true., !!!! & xmin,xmax, nopa, .false.) !!!! if (mnproc.eq.1) then !!!! write (nop,117) cname,nstep,time,0,coord,xmin,xmax !!!! call flush(nop) !!!! endif !1st tile !!!! cname = 'surty ' !!!! call zaiowr(surty,ip,.true., !!!! & xmin,xmax, nopa, .false.) !!!! if (mnproc.eq.1) then !!!! write (nop,117) cname,nstep,time,0,coord,xmin,xmax !!!! call flush(nop) !!!! endif !1st tile !!!! cname = 'sflice ' !!!! call zaiowr(util1,ip,.false., !mask on ice !!!! & xmin,xmax, nopa, .false.) !!!! if (mnproc.eq.1) then !!!! write (nop,117) cname,nstep,time,0,coord,xmin,xmax !!!! call flush(nop) !!!! endif !1st tile !!!! 117 format (a8,' =',i11,f11.3,i3,f7.3,1p2e16.7) !!!!c !!!! close (unit=nop) !!!! call zaiocl(nopa) !!!!c !!!! end subroutine Archive_ESMF !!!!#endif /* USE_ESMF */ subroutine HYCOM_Init #if defined (NERSC_VERSION) use mod_forcing_nersc use mod_random_forcing, only : randf, init_random_forcing, & rand_update use m_limits use m_icemodels_init #if defined (ICESTATE) use mod_icestate_init #endif #if defined (WAVES) ! use mod_restart ! use mod_common_wavesice use mod_waves_init, only: waves_setup #endif/* WAVES*/ #endif implicit none c c --- Initialize (before the 1st time step). c integer i,j,k,l,nm character*80 flnm,flnmra,flnmrb #if defined (NERSC_VERSION) character(len=80) tmparg integer*4, external :: iargc !type(year_info) :: rt #endif c include 'stmt_fns.h' c c c --- initialize SPMD processsing call xcspmd c c --- initialize timer names. c call xctmrn(40,'cnuity') call xctmrn(41,'tsadvc') call xctmrn(42,'momtum') call xctmrn(43,'barotp') call xctmrn(44,'thermf') call xctmrn(45,'ic****') call xctmrn(46,'mx****') call xctmrn(47,'conv**') call xctmrn(48,'diapf*') call xctmrn(49,'hybgen') call xctmrn(50,'trcupd') call xctmrn(51,'restrt') call xctmrn(52,'overtn') call xctmrn(53,'archiv') call xctmrn(54,'incupd') c c --- machine-specific initialization call machine c c --- initialize array i/o. call zaiost c c --- initiate named-pipe comparison utility call pipe_init c c --- initialize common variables c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! TODO (TW): could make these lines more efficient by going !! into blkdat, limits call blkdat #if defined (NERSC_VERSION) c --- Ensemble member from command-line arguments imem=1 if (iargc()==1) then C call getarg(1,tmparg) ; read(tmparg,'(i)') imem call getarg(1,tmparg) ; read(tmparg,*) imem ensflag=.true. elseif (iargc()>1) then call getarg(1,tmparg) ; read(tmparg,*) imem if (mnproc==1) then write(lp,*) 'too many arguments for model' write(lp,*) 'Usage: model ' end if call xcstop('(hycom)') stop '(hycom)' elseif (iargc()==0) then c --- New flag - means this is not an ensemble run ensflag=.false. end if if (mnproc==1) then write(lp,'(a)') '********************************' write(lp,'(a,i3)') 'Running HYCOM for member No ',imem write(lp,'(a)') '********************************' end if c --- Initialize ice variables in mod_common_ice call icedat c --- Initialize stuff specific for NERSC version call limits(yrflag, iniflg, relax, srelax, trelax, & tidflg,lbflag,nestfq,bnstfq,flxflg) #endif /* NERSC_VERSION */ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c !!!! nts_ice = icefrq !no. time steps between ice coupling !!!! nts_day = nint(86400.0d0/baclin) !no. time steps per day !!!! dsmall = baclin/86400.0d0 * 0.25d0 !1/4 of a time step in days !!!! dsmall2 = dsmall*2.0d0 !!!! if (dsurfq.ge.1.0) then !!!! ddsurf = dsurfq !!!! elseif (dsurfq.ne.0.0) then !!!! ddsurf = (baclin/86400.0d0)* !!!! & max(1,nint((86400.0d0*dsurfq)/baclin)) !!!! else !no surface archives !!!! ddsurf = 999999.0d0*(baclin/86400.0d0) !!!! endif !!!! if (diagfq.ge.1.0) then !!!! ddiagf = diagfq !!!! elseif (diagfq.ne.0.0) then !!!! ddiagf = (baclin/86400.0d0)* !!!! & max(1,nint((86400.0d0*diagfq)/baclin)) !!!! else !no 3-d archives !!!! ddiagf = huge*(baclin/86400.0d0) !!!! endif !!!! if (tilefq.ge.1.0) then !!!! dtilef = tilefq !!!! elseif (tilefq.ne.0.0) then !!!! dtilef = (baclin/86400.0d0)* !!!! & max(1,nint((86400.0d0*tilefq)/baclin)) !!!! else !no 3-d archives !!!! dtilef = huge*(baclin/86400.0d0) !!!! endif !!!! !!!! if (meanfq.ge.1.0) then !!!! dmeanf = meanfq !!!! elseif (meanfq.ne.0.0) then !!!! dmeanf = (baclin/86400.0d0)* !!!! & max(1,nint((86400.0d0*meanfq)/baclin)) !!!! else !no mean archives !!!! dmeanf = huge*(baclin/86400.0d0) !!!! endif !!!! if (rstrfq.eq.0.0) then ! no restart !!!! drstrf = rstrfq !!!! elseif (rstrfq.lt.0.0) then ! no restart at end of run !!!! drstrf = -rstrfq !!!! elseif (rstrfq.ge.1.0) then !!!! drstrf = rstrfq !!!! else !!!! drstrf = (baclin/86400.0d0)* !!!! & max(1,nint((86400.0d0*rstrfq)/baclin)) !!!! endif !!!! if (mnproc.eq.1) then !!!! write(lp,*) !!!! write(lp,*) 'ddsurf = ',ddsurf,nint((86400.0d0*ddsurf)/baclin) !!!! write(lp,*) 'ddiagf = ',ddiagf,nint((86400.0d0*ddiagf)/baclin) !!!! write(lp,*) 'dtilef = ',dtilef,nint((86400.0d0*dtilef)/baclin) !!!! write(lp,*) 'dmeanf = ',dmeanf,nint((86400.0d0*dmeanf)/baclin) !!!! write(lp,*) 'drstrf = ',drstrf,nint((86400.0d0*drstrf)/baclin) !!!! write(lp,*) !!!! write (lp,101) thkdf2,temdf2, !!!! & thkdf4, !!!! & veldf2,visco2, !!!! & veldf4,visco4, !!!! & diapyc,vertmx !!!! 101 format ( !!!! & ' turb. flux parameters:',1p/ !!!! & ' thkdf2,temdf2 =',2e9.2/ !!!! & ' thkdf4 =', e9.2/ !!!! & ' veldf2,visco2 =',2e9.2/ !!!! & ' veldf4,visco4 =',2e9.2/ !!!! & ' diapyc,vertmx =',2e9.2/) !!!! endif !1st tile c c --- days in year. c if (yrflag.eq.0) then c --- 360 days, starting Jan 16 dmonth = 30.0d0 dbimon = 60.0d0 dyear = 360.0d0 dyear0 = 0.0d0 elseif (yrflag.eq.1) then c --- 366 days, starting Jan 16 dmonth = 30.5d0 dbimon = 61.0d0 dyear = 366.0d0 dyear0 = 0.0d0 elseif (yrflag.eq.2) then c --- 366 days, starting Jan 1 c --- also implies high frequency atmospheric forcing dmonth = 30.5d0 dbimon = 61.0d0 dyear = 366.0d0 dyear0 = -15.0d0+dyear elseif (yrflag.eq.3) then c --- model day is calendar days since 01/01/1901 c --- also implies high frequency atmospheric forcing dyear = 365.25d0 dmonth = dyear/12.d0 dbimon = dyear/ 6.d0 dyear0 = -15.0d0+dyear else if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in hycom - unsupported yrflag value' write(lp,*) call flush(lp) endif !1st tile call xcstop('(hycom)') stop '(hycom)' !won't get here endif c c --- 'lstep' = number of barotropic time steps per baroclinic time step. c --- lstep m u s t be even. c !!!! lstep=nint(baclin/batrop) !!!! lstep=2*((lstep+1)/2) !!!! dlt=baclin/lstep !!!! if (mnproc.eq.1) then !!!! write (lp,'(i4,a/)') !!!! & lstep,' barotropic steps per baroclinic time step' !!!! endif !1st tile c c --- number of baroclinic time steps per day... !!!! nsteps_per_day = nint(86400.0/baclin) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!TW: need this!! c --- set up parameters defining the geographic environment c call geopar !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c #if defined (NERSC_VERSION) c --- Initialize random forcing !!!! if (randf) then !!!! call init_random_forcing() !!!! end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! TODO (TW): could make these lines more efficient by going !! into gridpinit, icemodels_init c --- Initialize gridp module call gridpinit() c c --- Initialize ice models call icemodels_init() c !!!!c --- Initialize tides (nersc version of tides) !!!! if (ltides) then !!!!c --- Time input is now hycom time, conversion is done in !!!!c --- tide_ast routine astr_time is set inside tides module !!!!c --- tide_ast init moved to first time step !!!! call tiderSB(onemm,1.0) !!!! if (mnproc==1) !!!! & write(lp,*) 'NERSC tides must be checked (new time input )' !!!! endif c #endif /*NERSC_VERSION*/ c c --- set up forcing functions c !!!! if (yrflag.lt.2) then !!!! call forfuna ! monthly atmospheric forcing !!!! endif !!!! if (jerlv0.eq.0) then !!!! call forfunk ! annual/monthly kpar !!!! endif !!!! call forfunp ! annual/monthly rivers !!!! !!!! call forfunr ! bimonthly/monthly climatology !!!! watcum=0. !!!! empcum=0. !!!!c !!!!c --- set minimum salinity for each isopycnic layer !!!! do k=2,kk !!!! cold=-3.0 !!!! salmin(k)=(sigma(k)-c1-cold*(c2+cold*(c4+c6*cold)))/ !!!! & (c3+cold*(c5+c7*cold)) !!!! enddo !!!!c !!!!c --- layer specific volume is defined as (1-theta)*thref !!!!c --- subtract constant 'thbase' from theta to reduce roundoff errors !!!!c !!!! if (vsigma) then !!!! call forfunv ! spacially varying isopycnal target densities !!!! else !!!! do k=1,kk !!!! theta(:,:,k)=sigma(k)-thbase !!!! enddo !!!! endif !!!!c !!!!c --- minimum depth of isopycnmal layers (pressure units). !!!!c !!!! if (isotop.lt.0.0) then !!!! call forfunt !spacially varying minimum depths !!!! else !!!! topiso(:,:)=onem*isotop !constant minimum depth !!!! endif !!!!c !!!!c --- tidal drag roughness (m/s) !!!!c !!!! if (drgscl.ne.0.0) then !!!! call forfund !tidal drag roughness !!!! else !!!! dragrh(:,:)=0.0 !!!! endif !!!!c !!!!c --- veldf2, veldf4 and thkdf4 may be spacially varying !!!!c !!!! call forfundf c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!TW: "limits stuff" - need this c --- model is to be integrated from time step 'nstep1' to 'nstep2' open( unit=uoff+99,file=trim(flnminp)//'limits') read( uoff+99,*) day1,day2 close(unit=uoff+99) c --- non-positive day1 indicates a new initialization, or c --- the start of a yrflag==3 case. linit =day1.le.0.0 day1 =abs(day1) if (mnproc==1) print*,'TW: day1,day2:',day1,day2 c dtime=day1 nstep1=nint(dtime*(86400.0d0/baclin)) dtime=nstep1/(86400.0d0/baclin) day1 =dtime c dtime=day2 nstep2=nint(dtime*(86400.0d0/baclin)) dtime=nstep2/(86400.0d0/baclin) day2 =dtime if (mnproc==1) print*,'TW2: day1,day2:',day1,day2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c !!!! if (mxlkpp) then !!!!c --- initialize kpp mixing !!!! call inikpp !!!! elseif (mxlmy) then !!!!c --- initialize m-y 2.5 mixing !!!! call inimy !!!! elseif (mxlgiss) then !!!!c --- initialize nasa giss mixing !!!! call inigiss !!!! endif c #if defined(USE_CCSM3) !!!! if (linit) then !!!! delt1 = baclin !!!! dhdx = c0p !!!! dhdy = c0p !!!! QICE = c0p !!!! AQICE = c0p !!!! QFLUX = c0p !!!!c !!!!c --- set up initial conditions !!!!c !!!!#if defined (NERSC_VERSION) !!!! call year_day(time,refyear,rt,yrflag) !!!!#endif !!!! mnth = imonth !!!! call inicon(mnth) !!!! trcrin = .false. !!!! call initrc(mnth) !!!!c !!!!c --- setup parameters defining tidal body forces !!!!c !!!! if (tidflg.gt.0) then !!!! time_8=dtime0 !'baroclinic' time for body force tides !!!! call tides_set(0) !!!! endif !!!! else !!!!c !!!!c --- start from restart file !!!!c !!!! delt1=baclin+baclin !!!! if (mnproc == 1) then !!!! open(1,file=trim(pointer_filename),form='formatted', !!!! & status='old') !!!! read(1,'(a)') flnmra !!!! read(1,'(a)') flnmrb !!!! read(1,*)time ! real, in days !!!! read(1,*)dtime ! real8, in days !!!! read(1,*)nstep,iyear,imonth,iday ! integer !!!! &, elapsed_days !!!! read(1,*)eom, eoy ! logical !!!! close(1) !!!! write(lp,*) ' flnm = ', flnmra !!!! write(lp,*) ' flnmh = ', flnmrb !!!! write(lp,*) ' pointer_filename = ', pointer_filename !!!! call flush (lp) !!!! write(lp,'(2a)')'(hycom) the restart info from: ', !!!! & trim(pointer_filename) !!!! write(lp,'(2a)')' restart file= ',trim(flnmra) !!!! write(lp,*) ' time= ',time !!!! write(lp,*) ' nstep= ',nstep !!!! write(lp,*) ' iyear= ',iyear !!!! write(lp,*) ' imonth= ',imonth !!!! write(lp,*) ' iday= ',iday !!!! write(lp,*) ' elapsed_days= ',elapsed_days !!!! write(lp,*) ' eom= ', eom !!!! write(lp,*) ' eoy= ', eoy !!!! call flush(lp) !!!! endif !mnproc==1 !!!! !broadcast the above info to all processors: !!!! call broadcast_scalar(flnmra, master_task) !!!! call broadcast_scalar(flnmrb, master_task) !!!! call broadcast_scalar(time, master_task) !!!! call broadcast_scalar(dtime, master_task) !!!! call broadcast_scalar(nstep, master_task) !!!! call broadcast_scalar(iyear, master_task) !!!! call broadcast_scalar(imonth, master_task) !!!! call broadcast_scalar(iday, master_task) !!!! call broadcast_scalar(elapsed_days, master_task) !!!! call broadcast_scalar(eom, master_task) !!!! call broadcast_scalar(eoy, master_task) !!!! !!!!!----------------------------------------------------------------------- !!!!! make any other needed adjustments after header file has been read !!!!! and report restart time in format recognized by ccsm3 !!!!!----------------------------------------------------------------------- !!!! call ccsm3_time_init1 !!!! !!!! call restart_in(nstep,dtime, flnmra, flnmrb) !!!! !!!! if (mnproc == 1 ) then !!!! write(lp,*)'(hycom) restartfile reading successful!' !!!! endif !!!! ! !!!! time0 = (iyear0-1)*days_in_year+days_in_prior_months(imonth0) !!!! & +iday0+0.0001 !!!! !the above makes time0 be the end of last integration (in days)! !!!! !!!! mnth = imonth !!!!#if defined (NERSC_VERSION) !!!! call year_day(time,refyear,rt,yrflag) !!!!#endif !!!! call initrc(mnth) !!!!c !!!!c --- setup parameters defining tidal body forces !!!!c !!!! if (tidflg.gt.0) then !!!! time_8=dtime0 !'baroclinic' time for body force tides !!!! call tides_set(0) !!!! endif !!!! endif !initial conditions #else if (linit .and. yrflag.lt.3) then !!!!c !!!!c --- set up initial conditions !!!!c !!!! nstep0=nstep1 !!!! dtime0=nstep0/(86400.0d0/baclin) !!!! time0=dtime0 !!!! delt1=baclin !!!! if (clmflg.eq.12) then !!!! mnth= 1.+nint(mod(dtime0+dyear0,dyear)/dmonth) !!!! elseif (clmflg.eq.6) then !!!! mnth=2*(1.+nint(mod(dtime0+dyear0,dyear)/dbimon))-1 !!!! endif !!!! call inicon(mnth) !!!! trcrin = .false. !!!!#if defined (NERSC_VERSION) !!!! call year_day(time,refyear,rt,yrflag) !!!!#endif !!!! call initrc(mnth) !!!!#if defined (NERSC_VERSION) !!!!c !!!!c --- This initializes the ice model, based on initial fields from inicon !!!! call iceinit !!!!#if defined (ICESTATE) !!!!c --- First attempt to initialize ICESTATE. Uses fields calc by iceinit !!!! call icestate_MICOM_init(hicem,ficem,hsnwm,tsrfm, !!!! & dp(:,:,1,1)/onem,saln(:,:,1,1),temp(:,:,1,1),rhosnw) !!!!#endif /* ICESTATE */ !!!!c !!!!#endif /* NERSC_VERSION */ !!!!c !!!!c --- setup parameters defining tidal body forces !!!!c !!!! if (tidflg.gt.0) then !!!! time_8=dtime0 !'baroclinic' time for body force tides !!!! call tides_set(0) !!!! endif !!!!c !!!!c --- output to archive file !!!!c !!!! m=mod(nstep0 ,2)+1 !!!! n=mod(nstep0+1,2)+1 !!!! nstep=nstep0 !!!! time=dtime0 !!!! call forday(dtime0,yrflag, iyear,jday,ihour) !!!!c !!!!!$OMP PARALLEL DO PRIVATE(j,k,l,i) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1,jj !!!! do l=1,isp(j) !!!! do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) !!!! tmix(i,j)=temp(i,j,1,n) !!!! smix(i,j)=saln(i,j,1,n) !!!! thmix(i,j)=th3d(i,j,1,n) !!!! enddo !i !!!! enddo !l !!!! if (isopyc .or. mxlkrt) then !!!! do l=1,isu(j) !!!! do i=max(1,ifu(j,l)),min(ii,ilu(j,l)) !!!! umix(i,j)=u(i,j,1,n) !!!! enddo !i !!!! enddo !l !!!! do l=1,isv(j) !!!! do i=max(1,ifv(j,l)),min(ii,ilv(j,l)) !!!! vmix(i,j)=v(i,j,1,n) !!!! enddo !i !!!! enddo !l !!!! endif !isopyc.or.mxlkrt !!!! enddo !j !!!!!$OMP END PARALLEL DO !!!!c !!!! if (mnproc.eq.1) then !!!! write (intvl,'(i3.3)') 0 !!!! endif !1st tile !!!! if (rstrfq.ne.0.0) then !don't write if benchmarking (no restart) !!!! call archiv(n, kk, iyear,jday,ihour, intvl) !!!! endif !!!!c else c c --- start from restart file c !!!!#if defined (NERSC_VERSION) #if defined (WAVES) !!initialise waves if (mnproc==1) then print *,'TW: in mod_hycom_wavesonly.F' print *,'day1=',day1 endif call waves_setup(day1)!!xcstop in here for WAVES_ONLY #endif/* WAVES */ !!!! if (mnproc==1) print*,'TW: calling read_restart_NERSC' !!!! call read_restart_NERSC(day1,nstep0,dtime0) !!!!Cc --- Initialize rt (date) !!!!C call year_day(day1,refyear,rt,yrflag) !!!!Cc --- Read one hycom restart file !!!!C call read_restart_mem(rt,nstep0,dtime0) !!!!Cc --- Read ice model(s) restart !!!!C call icemodels_read_restart(rt,nstep0,dtime0) !!!!Cc --- Read random restart !!!!C if (randf) then !!!!C call read_restart_rand_ab(rt) !!!!C end if !!!!#if defined (WEEKLY_AVERAGE) !!!!c --- Read weekly average restart !!!! call year_day(day1,refyear,rt,yrflag) !!!! call read_restart_ave(rt,baclin) !!!!#endif /* WEEKLY_AVERAGE */ c !!!!#else /* ! NERSC_VERSION */ !!!! flnmra = trim(flnmrsi)//'.a' !!!! flnmrb = trim(flnmrsi)//'.b' !!!! call restart_in(nstep0,dtime0, flnmra,flnmrb) !!!!#endif /* ! NERSC_VERSION */ c !!!! if (linit) then !!!!c !!!!c --- start a new calendar-day (yrflag==3) case. !!!! if (mnproc.eq.1) then !!!! time0=dtime0 !!!! write (lp,'(9x,a,f8.1,a,f8.1,a,i9)') !!!! & 'restart file for day',time0, !!!! & ', is treated as day', day1,' step',nstep1 !!!! endif !1st tile !!!! nstep0=nstep1 !!!! dtime0=nstep0/(86400.0d0/baclin) !!!! time0=dtime0 !!!! delt1=baclin+baclin !!!! else !!!! nstep0=nint(dtime0*(86400.0d0/baclin)) !!!! dtime0=nstep0/(86400.0d0/baclin) !!!! time0=dtime0 !!!! delt1=baclin+baclin !!!! if (mnproc.eq.1) then !!!! write (lp,'(a,f8.1,a,i9,a, a,f8.1,a,i9,a)') !!!! & 'restart on day',time0,' (step',nstep0,')', !!!! & ', wanted day', day1,' (step',nstep1,')' !!!! endif !1st tile !!!! if (nstep0.ne.nstep1) then !!!! if (mnproc.eq.1) then !!!! write(lp,'(/a/a,f8.1/a,f8.1/)') !!!! & 'error in hycom - wrong restart (or limits) file', !!!! & 'restart file day is ',time0, !!!! & 'limits start day is ',day1 !!!! endif !1st tile !!!! call xcstop('(hycom)') !!!! stop '(hycom)' !won't get here !!!! endif !nstep0.ne.nstep1 !!!! endif !linit:else !!!!c !!!! if (clmflg.eq.12) then !!!! mnth= 1.+nint(mod(dtime0+dyear0,dyear)/dmonth) !!!! elseif (clmflg.eq.6) then !!!! mnth=2*(1.+nint(mod(dtime0+dyear0,dyear)/dbimon))-1 !!!! endif !!!! call initrc(mnth) !!!!c !!!!c --- setup parameters defining tidal body forces !!!!c !!!! if (tidflg.gt.0) then !!!! time_8=dtime0 !'baroclinic' time for body force tides !!!! call tides_set(0) !!!! endif !!!!c !!!! if (trcout .and. .not.trcrin) then !!!!c !!!!c --- new tracers, so output to archive file !!!!c !!!! m=mod(nstep0 ,2)+1 !!!! n=mod(nstep0+1,2)+1 !!!! call momtum_hs(n,m) !calculate srfhgt !!!! nstep=nstep0 !!!! time=dtime0 !!!! call forday(dtime0,yrflag, iyear,jday,ihour) !!!! if (mnproc.eq.1) then !!!! write (intvl,'(i3.3)') 0 !!!! endif !1st tile !!!! call archiv(n, kk, iyear,jday,ihour, intvl) !!!! endif !archive output endif !initial conditions #endif /* USE_CCSM3:else */ c c --- set barotp.pot.vort. and layer thickness (incl.bottom pressure) at c --- u,v points c !!!! call dpthuv !!!!c !!!! call xctilr(dp( 1-nbdy,1-nbdy,1,1),1,2*kk, nbdy,nbdy, halo_ps) !!!! call xctilr(dpmixl(1-nbdy,1-nbdy, 1),1, 2, nbdy,nbdy, halo_ps) !!!! call xctilr(thkk( 1-nbdy,1-nbdy,1 ),1, 2, nbdy,nbdy, halo_ps) !!!! call xctilr(psikk( 1-nbdy,1-nbdy,1 ),1, 2, nbdy,nbdy, halo_ps) !!!!c !!!! margin = nbdy !!!!c !!!! do nm=1,2 !!!!c !!!!!$OMP PARALLEL DO PRIVATE(j,k,l,i) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1-margin,jj+margin !!!! do l=1,isp(j) !!!! if (nm.eq.mod(nstep+1,2)+1) then !!!! do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) !!!! dpbl(i,j)=dpmixl(i,j,nm) !!!! enddo !!!! endif !!!! do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) !!!! p(i,j,1)=0.0 !!!! do k=1,kk !!!! p(i,j,k+1)=p(i,j,k)+dp(i,j,k,nm) !!!! enddo !!!! enddo !!!! enddo !!!! enddo !!!!!$OMP END PARALLEL DO !!!!c !!!! call dpudpv(dpu(1-nbdy,1-nbdy,1,nm), !!!! & dpv(1-nbdy,1-nbdy,1,nm), !!!! & p,depthu,depthv, max(0,margin-1)) !!!!c !!!! if (.false.) then !!!!c !!!!c --- ISOPYC TO HYBRID RESTART ONLY !!!! nstep=nstep1 !!!! n=nm !!!! m=mod(n,2)+1 !!!! call hybgen(m,n) !!!! call xctilr(dp(1-nbdy,1-nbdy,1,n),1,kk, nbdy,nbdy, halo_ps) !!!! margin = nbdy !!!!c !!!!!$OMP PARALLEL DO PRIVATE(j,k,l,i) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1-margin,jj+margin !!!! do l=1,isp(j) !!!!!DIR$ PREFERVECTOR !!!! do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) !!!! p(i,j,1)=0.0 !!!! do k=1,kk !!!! p(i,j,k+1)=p(i,j,k)+dp(i,j,k,n) !!!! enddo !!!! enddo !!!! enddo !!!! enddo !!!!!$OMP END PARALLEL DO !!!! call dpudpv(dpu(1-nbdy,1-nbdy,1,n), !!!! & dpv(1-nbdy,1-nbdy,1,n), !!!! & p,depthu,depthv, max(0,margin-1)) !!!! call pipe_comparall(m,n, 'hybgen, step') !!!! endif !isopyc to hybrid restart only !!!!c !!!! enddo !nm=1,2 c c --- mean archive initialization. c !!!! if (meanfq.ne.0.0) then !!!! call mean_allocate !!!! m=mod(nstep0 ,2)+1 !!!! n=mod(nstep0+1,2)+1 !!!! nstep=nstep0 !!!! time =dtime0 !!!! call mean_zero(time) !!!! call momtum_hs(n,m) !calculate srfhgt !!!! if (tidflg.gt.0 .and. meanfq.eq.1.0) then !!!! call mean_add(n, 1.0) !25-hour average daily !!!! else !!!! call mean_add(n, 0.5) !!!! endif !!!! endif c !!!!#if ! defined(USE_CCSM3) nod=14 nstep=nstep1 if (mnproc.eq.1) then write(lp ,'(/2(a,f8.1),2(a,i9),a/)') 'model starts at day', & day1,', goes to day',day2 open (unit=nod,file='summary_out',status='unknown') write(nod,'(/2(a,f8.1),2(a,i9),a/)') 'model starts at day', & day1,', goes to day',day2 endif !1st tile !!!!#endif !!!!c !!!! timav=time0 !!!! m=mod(nstep ,2)+1 !!!! n=mod(nstep+1,2)+1 !!!!c !!!! call pipe_comparall(m,n, 'restrt, step') !!!!c !!!! if (synflt) then !!!!c --- initialize synthetic floats/moorings !!!! call floats_init(m,n,time0) !!!! margin = nbdy !!!! endif !!!!#if defined(USE_CCSM3) !!!!!---------------------------------------------------------- !!!!! in the coupled model, forcing is received from the !!!!! coupler; therefore code for reading stand-alone version !!!!! forcing fields has been removed !!!!!---------------------------------------------------------- !!!! !!!! nstep0=nstep !!!! m=mod(nstep ,2)+1 !!!! n=mod(nstep+1,2)+1 !!!! !!!! k1m=m !NOT k1m=1+mm !!!! k1n=n !NOT k1n=1+nn !!!! !!!! !!!! !------------------------------------------------------------- !!!! ! initialize namelists, grid information, and initial message- !!!! ! passing for cpl6 !!!! !------------------------------------------------------------- !!!! call ccsm3_init_coupled(flnmgrdd) !!!! !!!! if (mnproc == 1) then !!!! write(lp,*)'(hycom) ccsm3_init_coupled completed' !!!! call flush(lp) !!!! endif !!!! !!!! !!!! if (mnproc == 1) then !!!! write (*,*) '(hycom) model starts from day, nstep=', iday,nstep !!!! endif !!!! !!!! l0=1 !!!! l1=2 !!!! l2=3 !!!! l3=4 !!!! w0=0.0 !!!! w1=1.0 !!!! w2=0.0 !!!! w3=0.0 !!!!c !!!! if (jerlv0.eq.0) then !!!!c --- read in kpar field for 4 consecutive months !!!! mk1=imonth !!!! mk0=mod(mk1+10,12)+1 !!!! mk2=mod(mk1, 12)+1 !!!! mk3=mod(mk2, 12)+1 !!!! lk0=1 !!!! lk1=2 !!!! lk2=3 !!!! lk3=4 !!!! call rdkpar(mk0,lk0) !!!! call rdkpar(mk1,lk1) !!!! call rdkpar(mk2,lk2) !!!! call rdkpar(mk3,lk3) !!!! endif !!!!c !!!! if (clmflg.eq.12) then !!!!c --- read in relaxation climatology fields for 4 consecutive months !!!! mc1=imonth !!!! mc0=mod(mc1+10,12)+1 !!!! mc2=mod(mc1, 12)+1 !!!! mc3=mod(mc2, 12)+1 !!!! lc0=1 !!!! lc1=2 !!!! lc2=3 !!!! lc3=4 !!!! call rdrlax(mc0,lc0) !!!! call rdrlax(mc1,lc1) !!!! call rdrlax(mc2,lc2) !!!! call rdrlax(mc3,lc3) !!!! elseif (clmflg.eq.6) then !!!!c --- read in relaxation fields for 4 consecutive bi-months !!!! mc1=imonth !!!! mc0=mod(mc1+4,6)+1 !!!! mc2=mod(mc1, 6)+1 !!!! mc3=mod(mc2, 6)+1 !!!! lc0=1 !!!! lc1=2 !!!! lc2=3 !!!! lc3=4 !!!! call rdrlax(2*mc0-1,lc0) !!!! call rdrlax(2*mc1-1,lc1) !!!! call rdrlax(2*mc2-1,lc2) !!!! call rdrlax(2*mc3-1,lc3) !!!! else !!!! if (mnproc.eq.1) then !!!! write(lp,'(/ a /)') 'error in hycom - unsupported clmflg value' !!!! call flush(lp) !!!! endif !1st tile !!!! call xcstop('(hycom)') !!!! stop '(hycom)' !won't get here !!!! endif !!!! !!!! nod = 14 !!!! write(cyear,'(i4.4)') iyear !!!!Cpg - Write out summary files for each month !!!! write(flnm,'(a,a,i4.4,a,i2.2)') !!!! & trim(flnmsumd),'summary_out.',iyear,'-',imonth !!!! if (mnproc==1) then !!!! open(unit=nod,file=flnm,status='unknown') !!!! write(nod,'(3a,i2.2)') !!!! & 'Summary of diagnostic outputs for year ',cyear !!!! & , ' month ', imonth !!!!Cpg !!!! write(nod,*) !!!! call flush(nod) !!!! endif !!!!#else !!!!c !!!! if (yrflag.lt.2) then !!!!c !!!!c --- read in forcing fields for 4 consecutive months !!!! ma1=1.+mod(dtime0+dyear0,dyear)/dmonth !!!! ma0=mod(ma1+10,12)+1 !!!! ma2=mod(ma1, 12)+1 !!!! ma3=mod(ma2, 12)+1 !!!! l0=1 !!!! l1=2 !!!! l2=3 !!!! l3=4 !!!! call rdforf(ma0,l0) !!!! call rdforf(ma1,l1) !!!! call rdforf(ma2,l2) !!!! call rdforf(ma3,l3) !!!! else !!!!c !!!!c --- initial day of high frequency atmospheric forcing. !!!!c --- only two fields are used (linear interpolation in time). !!!! l0=1 !!!! l1=2 !!!! l2=3 !!!! l3=4 !!!! if (windf) then !!!! w0=-99.9 !!!! w1=-99.0 !!!! w2=0.0 !!!! w3=0.0 !!!!#if defined (NERSC_VERSION) && defined (FORCING_INLINE) !!!! call nersc_forcing_inline( !!!! & dtime0,1.+mod(dtime0+dyear0,dyear)/dmonth) !!!!#else !!!! call forfunh(dtime0) !!!!#endif !!!!#if defined (NERSC_VERSION) !!!! if (randf) then !!!! call rand_update(dtime0,1,.true.) !!!! call rand_update(dtime0,2,.true.) !!!! end if !!!!#endif !!!! else !!!! w0=0.0 !!!! w1=0.0 !!!! w2=0.0 !!!! w3=0.0 !!!! endif !!!! endif !!!!c !!!! if (jerlv0.eq.0) then !!!!c --- read in kpar field for 4 consecutive months !!!! mk1=1.+mod(dtime0+dyear0,dyear)/dmonth !!!! mk0=mod(mk1+10,12)+1 !!!! mk2=mod(mk1, 12)+1 !!!! mk3=mod(mk2, 12)+1 !!!! lk0=1 !!!! lk1=2 !!!! lk2=3 !!!! lk3=4 !!!! call rdkpar(mk0,lk0) !!!! call rdkpar(mk1,lk1) !!!! call rdkpar(mk2,lk2) !!!! call rdkpar(mk3,lk3) !!!! endif !!!!c !!!! if (priver) then !!!!c --- read in rivers field for 4 consecutive months !!!! mr1=1.+mod(dtime0+dyear0,dyear)/dmonth !!!! mr0=mod(mr1+10,12)+1 !!!! mr2=mod(mr1, 12)+1 !!!! mr3=mod(mr2, 12)+1 !!!! lr0=1 !!!! lr1=2 !!!! lr2=3 !!!! lr3=4 !!!! call rdrivr(mr0,lr0) !!!! call rdrivr(mr1,lr1) !!!! call rdrivr(mr2,lr2) !!!! call rdrivr(mr3,lr3) !!!! endif !!!!c !!!! if (clmflg.eq.12) then !!!!c --- read in relaxation climatology fields for 4 consecutive months !!!! mc1=1.+mod(dtime0+dyear0,dyear)/dmonth !!!! mc0=mod(mc1+10,12)+1 !!!! mc2=mod(mc1, 12)+1 !!!! mc3=mod(mc2, 12)+1 !!!! lc0=1 !!!! lc1=2 !!!! lc2=3 !!!! lc3=4 !!!! call rdrlax(mc0,lc0) !!!! call rdrlax(mc1,lc1) !!!! call rdrlax(mc2,lc2) !!!! call rdrlax(mc3,lc3) !!!! elseif (clmflg.eq.6) then !!!!c --- read in relaxation fields for 4 consecutive bi-months !!!! mc1=1.+mod(dtime0+dyear0,dyear)/dbimon !!!! mc0=mod(mc1+4,6)+1 !!!! mc2=mod(mc1, 6)+1 !!!! mc3=mod(mc2, 6)+1 !!!! lc0=1 !!!! lc1=2 !!!! lc2=3 !!!! lc3=4 !!!! call rdrlax(2*mc0-1,lc0) !!!! call rdrlax(2*mc1-1,lc1) !!!! call rdrlax(2*mc2-1,lc2) !!!! call rdrlax(2*mc3-1,lc3) !!!! else !!!! if (mnproc.eq.1) then !!!! write(lp,'(/ a /)') 'error in hycom - unsupported clmflg value' !!!! call flush(lp) !!!! endif !1st tile !!!! call xcstop('(hycom)') !!!! stop '(hycom)' !won't get here !!!! endif !!!!c !!!! if (bnstfq.ne.0.0) then ! initialize barotropic boundary input !!!! wb0=-99.0 !!!! wb1=-99.0 !!!! call rdbaro(dtime0) !!!! endif !!!!c !!!! if (nestfq.ne.0.0) then ! initialise 3-d nesting input !!!! wn0=-99.0 !!!! wn1=-99.0 !!!! call rdnest(dtime0) !!!! endif !!!!#endif /* USE_CCSM3:else */ c c --- initialize incremental update. c !!!! if (incflg.ne.0) then !!!! call incupd_init(dtime0) !!!!c !!!! if (incstp.eq.1) then !!!! call xctmr0(54) !!!! call incupd(1) !!!! call incupd(2) !!!! call xctmr1(54) !!!! endif ! full insertion of update !!!! endif !!!!#if defined(USE_ESMF) !!!!c !!!!c --- Fill the Export State with the initial fields. !!!!c !!!! m=mod(nstep0 ,2)+1 !!!! n=mod(nstep0+1,2)+1 !!!! call momtum_hs(n,m) !!!! call Export_ESMF !!!!#else !!!!c !!!!c --- Only here for compatibility with coupled runs. !!!!c !!!! m=mod(nstep0 ,2)+1 !!!! n=mod(nstep0+1,2)+1 !!!! call momtum_hs(n,m) !!!!#endif !!!!c !!!!c --- report initialization time. !!!!c !!!! call xctmrp c end subroutine HYCOM_Init subroutine HYCOM_Run #if defined(NERSC_VERSION) use mod_forcing_waves, only: run_waves_only !!!! use mod_nesting, only : lnesto, lnesti, !!!! & nestcondsave, nestdto, nestdti !!!! use mod_random_forcing, only : randf, rand_update !!!!C use m_icemodels_save_restart !!!! use m_icemodels_step !!!!#if defined (DAILY_AVERAGE) !!!! use mod_daily_average !!!!#endif #endif /* NERSC_VERSION */ implicit none !!!!#if defined(USE_ESMF) !!!!c !!!!c --- Calling parameters !!!! type(ESMF_GridComp) :: gridComp !!!! type(ESMF_State) :: impState !!!! type(ESMF_State) :: expState !!!! type(ESMF_Clock) :: extClock !!!! integer :: rc !!!!#endif c c --- ------------------------- c --- execute a single timestep c --- ------------------------- c !!!! logical lfatal !!!! integer i,j,k,ktr,l,nm !!!! character*80 flnmra,flnmrb !!!!c !!!! logical hycom_isnaninf !function to detect NaN and Inf !!!!c !!!!#if defined (NERSC_VERSION) !!!! integer :: ierr !!!! !integer :: outday, ierr !!!! !logical :: loutput !!!!#endif c !!!! include 'stmt_fns.h' c c --- letter 'm' refers to mid-time level (example: dp(i,j,k,m) ) c --- letter 'n' refers to old and new time level c if (mnproc==1) print*, & 'TW (mod_hycom_wavesonly.F): running waves-only' call run_waves_only end_of_run = .true. if (mnproc==1) print*, & 'TW (mod_hycom_wavesonly.F): waves-only run finished' !!call xcstop('TW (mod_hycom_wavesonly.F): waves-only run finished') !!stop !!!! m=mod(nstep ,2)+1 !!!! n=mod(nstep+1,2)+1 !!!!c !!!!#if defined(USE_CCSM3) !!!! k1m=m !NOT k1m=1+mm (DBI) !!!! k1n=n !NOT k1n=1+nn (DBI) !!!! !!!!!------------------------------------------------------------------------- !!!!! exchange fields and fluxes with the ccsm3 flux coupler !!!!!------------------------------------------------------------------------- !!!! call ccsm3_set_coupled_forcing !!!! !!!!!-------------------------------------------------------- !!!!! stop if cpl has sent the stop signal !!!!!-------------------------------------------------------- !!!! end_of_run = stop_now == 1 !!!! if (end_of_run) then !!!!c --- output float restart file !!!! if (synflt) then !!!! call floats_restart !!!! endif !synflt !!!! return !!!! endif !!!! !!!!!------------------------------------------------------------------------- !!!!! ccsm3-related time-keeping !!!!!------------------------------------------------------------------------- !!!! call ccsm3_time_advance (time, dtime) !!!!c !!!! hisurf=.false. !!!! histry=.false. !!!! hitile=.false. !!!! histmn=dohist !!!! restrt=dorestart .or. (cpl_write_restart .and. eod) !!!! diagno=mod(dtime+dsmall,ddiagf).lt.dsmall2 .or. restrt !!!!c !!!!c --- set weights for quasi-hermite time interpolation for kpar. !!!! if (jerlv0.eq.0) then !!!!c --- monthly fields. !!!! if (imonth.ne.mk1) then !!!! mk1=imonth !!!! mk0=mod(mk1+10,12)+1 !!!! mk2=mod(mk1, 12)+1 !!!! mk3=mod(mk2, 12)+1 !!!! lt =lk0 !!!! lk0=lk1 !!!! lk1=lk2 !!!! lk2=lk3 !!!! lk3=lt !!!! call rdkpar(mk3,lk3) !!!! endif !!!! wk1=1.0 !constant for the month !!!! wk2=0.0 !!!! wk0=0.0 !!!! wk3=0.0 !!!! endif !!!!c !!!!c --- set weights for quasi-hermite time interpolation for temperature, !!!!c --- salinity and pressure relaxation fields. !!!! if (clmflg.eq.12) then !!!!c --- monthly fields. !!!! if (imonth.ne.mc1) then !!!! mc1=imonth !!!! mc0=mod(mc1+10,12)+1 !!!! mc2=mod(mc1, 12)+1 !!!! mc3=mod(mc2, 12)+1 !!!! lt =lc0 !!!! lc0=lc1 !!!! lc1=lc2 !!!! lc2=lc3 !!!! lc3=lt !!!! call rdrlax(mc3,lc3) !!!! endif !!!! wc1=1.0 !constant for the month !!!! wc2=0.0 !!!! wc0=0.0 !!!! wc3=0.0 !!!! endif !!!!cdiag if (mnproc.eq.1) then !!!!cdiag write (lp,'(i9,'' relax time interpolation: months'',4i3, !!!!cdiag. '', weights '',4f6.3)') nstep,lc0,lc1,lc2,lc3,wc0,wc1,wc2,wc3 !!!!cdiag endif !1st tile !!!!#else !!!! nstep=nstep+1 !!!! dtime=dtime0+(nstep-nstep0)/(86400.0d0/baclin) !!!!c !!!! time =dtime !!!! time_8=dtime !'baroclinic' time for body force tides !!!! if (tidflg.gt.0 .and. !!!! & mod(dtime+dsmall,hours1).lt.dsmall2) then !!!! call tides_detide(n, .true.) !update 25-hour average !!!! endif !!!! hisurf=mod(dtime+dsmall,ddsurf).lt.dsmall2 !!!! histry=mod(dtime+dsmall,ddiagf).lt.dsmall2 !!!! hitile=mod(dtime+dsmall,dtilef).lt.dsmall2 !!!! histmn=mod(dtime+dsmall,dmeanf).lt.dsmall2 .or. !!!! & nstep.ge.nstep2 !!!! if (rstrfq.eq.0.0) then ! no restart !!!! restrt=.false. ! for benchmark cases only !!!! elseif (rstrfq.lt.0.0) then ! no restart at end of run !!!! restrt=mod(dtime-dtime0+dsmall,drstrf).lt.dsmall2 !!!! else !!!! restrt=mod(dtime-dtime0+dsmall,drstrf).lt.dsmall2 .or. !!!! & nstep.ge.nstep2 !!!! endif !!!! diagno=mod(dtime+dsmall,ddiagf).lt.dsmall2 .or. !!!! & restrt .or. nstep.ge.nstep2 !!!! if (yrflag.lt.2) then !!!!c !!!!c --- set weights for quasi-hermite time interpolation for !!!!c --- monthly atmospheric forcing fields !!!! x=1.+mod(dtime+dyear0,dyear)/dmonth !!!!c --- keep quadruplet of forcing functions centered on model time !!!! if (int(x).ne.ma1) then !!!! ma1=x !!!! ma0=mod(ma1+10,12)+1 !!!! ma2=mod(ma1, 12)+1 !!!! ma3=mod(ma2, 12)+1 !!!! lt=l0 !!!! l0=l1 !!!! l1=l2 !!!! l2=l3 !!!! l3=lt !!!!c --- newest set of forcing functions overwrites set no longer needed !!!! call rdforf(ma3,l3) !!!! endif !!!! x=mod(x,1.) !!!! x1=1.-x !!!! w1=x1*(1.+x *(1.-1.5*x )) !!!! w2=x *(1.+x1*(1.-1.5*x1)) !!!! w0=-.5*x *x1*x1 !!!! w3=-.5*x1*x *x !!!!cdiag if (mnproc.eq.1) then !!!!cdiag write (lp,'(i9,'' atmos time interpolation: months'',4i3, !!!!cdiag. '', weights '',4f6.3)') nstep,l0,l1,l2,l3,w0,w1,w2,w3 !!!!cdiag endif !1st tile !!!! elseif (windf) then !!!!c !!!!c --- set weights and fields for high frequency atmospheric forcing. !!!!c --- only two fields are used (linear interpolation in time). !!!!#if defined (NERSC_VERSION) && defined (FORCING_INLINE) !!!! call nersc_forcing_inline( !!!! & dtime, 1.+mod(dtime+dyear0,dyear)/dmonth) !!!!#else !!!! call forfunh(dtime) !!!!#endif !!!!#if defined (NERSC_VERSION) !!!! if (randf) call rand_update(dtime ,2,.false.) !!!!#endif !!!! endif !!!!c !!!!c --- set weights for quasi-hermite time interpolation for kpar. !!!! if (jerlv0.eq.0) then !!!!c --- monthly fields. !!!! x=1.+mod(dtime+dyear0,dyear)/dmonth !!!! if (int(x).ne.mk1) then !!!! mk1=x !!!! mk0=mod(mk1+10,12)+1 !!!! mk2=mod(mk1, 12)+1 !!!! mk3=mod(mk2, 12)+1 !!!! lt =lk0 !!!! lk0=lk1 !!!! lk1=lk2 !!!! lk2=lk3 !!!! lk3=lt !!!! call rdkpar(mk3,lk3) !!!! endif !!!! x=mod(x,1.) !!!! x1=1.-x !!!! wk1=x1*(1.+x *(1.-1.5*x )) !!!! wk2=x *(1.+x1*(1.-1.5*x1)) !!!! wk0=-.5*x *x1*x1 !!!! wk3=-.5*x1*x *x !!!! endif !!!!c !!!!c --- set weights for quasi-hermite time interpolation for rivers. !!!! if (priver) then !!!!c --- monthly fields. !!!! x=1.+mod(dtime+dyear0,dyear)/dmonth !!!! if (int(x).ne.mr1) then !!!! mr1=x !!!! mr0=mod(mr1+10,12)+1 !!!! mr2=mod(mr1, 12)+1 !!!! mr3=mod(mr2, 12)+1 !!!! lt =lr0 !!!! lr0=lr1 !!!! lr1=lr2 !!!! lr2=lr3 !!!! lr3=lt !!!! call rdrivr(mr3,lr3) !!!! endif !!!! x=mod(x,1.) !!!! x1=1.-x !!!! wr1=x1*(1.+x *(1.-1.5*x )) !!!! wr2=x *(1.+x1*(1.-1.5*x1)) !!!! wr0=-.5*x *x1*x1 !!!! wr3=-.5*x1*x *x !!!! endif !!!!c !!!!c --- set weights for quasi-hermite time interpolation for temperature, !!!!c --- salinity and pressure relaxation fields. !!!! if (clmflg.eq.12) then !!!!c --- monthly fields. !!!! x=1.+mod(dtime+dyear0,dyear)/dmonth !!!! if (int(x).ne.mc1) then !!!! mc1=x !!!! mc0=mod(mc1+10,12)+1 !!!! mc2=mod(mc1, 12)+1 !!!! mc3=mod(mc2, 12)+1 !!!! lt =lc0 !!!! lc0=lc1 !!!! lc1=lc2 !!!! lc2=lc3 !!!! lc3=lt !!!! call rdrlax(mc3,lc3) !!!! endif !!!! x=mod(x,1.) !!!! x1=1.-x !!!! wc1=x1*(1.+x *(1.-1.5*x )) !!!! wc2=x *(1.+x1*(1.-1.5*x1)) !!!! wc0=-.5*x *x1*x1 !!!! wc3=-.5*x1*x *x !!!! elseif (clmflg.eq.6) then !!!!c --- bi-monthly fields. !!!! x=1.+mod(dtime+dyear0,dyear)/dbimon !!!! if (int(x).ne.mc1) then !!!! mc1=x !!!! mc0=mod(mc1+4,6)+1 !!!! mc2=mod(mc1, 6)+1 !!!! mc3=mod(mc2, 6)+1 !!!! lt =lc0 !!!! lc0=lc1 !!!! lc1=lc2 !!!! lc2=lc3 !!!! lc3=lt !!!! call rdrlax(2*mc3-1,lc3) !!!! endif !!!! x=mod(x,1.) !!!! x1=1.-x !!!! wc1=x1*(1.+x *(1.-1.5*x )) !!!! wc2=x *(1.+x1*(1.-1.5*x1)) !!!! wc0=-.5*x *x1*x1 !!!! wc3=-.5*x1*x *x !!!! endif !!!!cdiag if (mnproc.eq.1) then !!!!cdiag write (lp,'(i9,'' relax time interpolation: months'',4i3, !!!!cdiag. '', weights '',4f6.3)') nstep,lc0,lc1,lc2,lc3,wc0,wc1,wc2,wc3 !!!!cdiag endif !1st tile !!!!#endif /* USE_CCSM3:else */ !!!!c !!!!#if defined(USE_ESMF) !!!! if (get_import) then ! new ESMF Import fields !!!! call Import_ESMF !!!! if (nstep.eq.nstep0+1) then !!!! time=dtime0 !!!! call forday(dtime0,yrflag, iyear,jday,ihour) !!!! call Archive_ESMF(iyear,jday,ihour) !!!! time=dtime !!!! endif !initial ESMF Archive !!!! endif !get_import !!!!#endif !!!!c !!!! if (bnstfq.ne.0.0) then ! new fields for baro nesting !!!! call rdbaro(dtime) !!!! endif !!!!c !!!! if (nestfq.ne.0.0) then ! new fields for 3-d nesting !!!! call rdnest(dtime) !!!! endif !!!!#if defined (NERSC_VERSION) !!!! if (ltides) then !!!! if (nstep==nstep1+1) then !!!! call tide_ast(dtime,yrflag) !!!! elseif (dtime >= (astr_time+0.5*astr_dt)) then !!!! call tide_ast(astr_time+astr_dt,yrflag) !!!! endif !!!! if (mnproc==1) !!!! & write(lp,*) 'NERSC tides must be checked (new time input )' !!!! nersctide_time=dtime ! Set here, used in latbdf !!!! endif !!!!#endif /* NERSC_VERSION */ !!!!c !!!! call pipe_comparall(m,n, 'ENTER , step') !!!! call xctmr0(40) !!!! call cnuity(m,n) !!!! call xctmr1(40) !!!! call pipe_comparall(m,n, 'cnuity, step') !!!! call xctmr0(41) !!!! call tsadvc(m,n) !!!! call xctmr1(41) !!!! call pipe_comparall(m,n, 'tsadvc, step') !!!! call xctmr0(42) !!!! if (momtyp.eq.2) then !!!! call momtum(m,n) !!!! else !momtyp.eq.4 !!!! call momtum4(m,n) !!!! endif !!!! call xctmr1(42) !!!! call pipe_comparall(m,n, 'momtum, step') !!!! call xctmr0(43) !!!! call barotp(m,n) !!!! call xctmr1(43) !!!! call pipe_comparall(m,n, 'barotp, step') !!!! call xctmr0(44) !!!!#if defined(USE_CCSM3) !!!! call thermf_c(m,n) !,mm,nn,k1m,k1n) !!!! call xctmr1(44) !!!! call pipe_comparall(m,n, 'thermf, step') !!!! !get potential freezing/melting heat flux needed for ice model. !!!! !mixed layer T and S are adjusted only TWICE in each coupling interval. !!!! ! !!!! call ice_formation(k1n, ice_ts) !!!! !=================================================================== !!!! !if (.false. .and. ice_ts) then !DBI: switched this call off! !!!! ! call ice_flx_to_coupler ! work done in ice_formation now! !!!! !endif !!!!#else !!!!#if defined(NERSC_VERSION) !!!! if (flxflg.eq.99) then !!!! call year_day(time,refyear,rt,yrflag) !!!! call icemodels_step(m,n,nstep==nstep1+1,rt,rforce,dtime) !!!!c !!!!c --- Note that thermf takes care of relaxation / baroclinic nesting !!!!c --- NERSC_VERSION cpp statements are included in that file !!!! else !!!! if (mnproc==1) !!!! & write(lp,'(a)') 'Warning: NERSC_VERSION, but flxflg.ne.99' !!!! endif !!!!#endif !!!! call thermf(m,n) !!!! call xctmr1(44) !!!! call pipe_comparall(m,n, 'thermf, step') !!!! if (icegln) then !!!! call xctmr0(45) !!!! call icloan(m,n) !!!! call thermf_oi(m,n) !!!! call xctmr1(45) !!!! call pipe_comparall(m,n, 'icloan, step') !!!! elseif (iceflg.ge.2) then !!!! call xctmr0(45) !!!! call thermf_oi(m,n) !!!! call xctmr1(45) !!!! call pipe_comparall(m,n, 'icecpl, step') !!!! else !!!! call xctmr0(45) !!!! call thermf_oi(m,n) !!!! call xctmr1(45) !!!! call pipe_comparall(m,n, 'thermf, step') !thermf_oi !!!! endif !icegln:icecpl:else !!!!#endif /* USE_CCSM3:else */ !!!! if (trcout) then !!!! call xctmr0(50) !!!! call trcupd(m,n) !!!! call xctmr1(50) !!!! call pipe_comparall(m,n, 'trcupd, step') !!!! endif !trcout !!!! if (hybrid) then !!!! diagsv = diagno !!!! diagno = diagno .or. nstep.eq.nstep0+1 .or. !!!! & histry .or. hisurf .or. hitile .or. !!!! & mod(dtime+dsmall,days6).lt.dsmall2 !!!! if (mxlkpp .or. mxlmy .or. mxlgiss) then !!!! call xctmr0(46) !!!! call mxkprf(m,n) !!!! call xctmr1(46) !!!! call pipe_comparall(m,n, 'mxkprf, step') !!!! elseif (mxlpwp) then !!!! call xctmr0(46) !!!! call mxpwp(m,n) !!!! call xctmr1(46) !!!! call pipe_comparall(m,n, 'mxpwp, step') !!!! elseif (mxlkta) then !!!! call xctmr0(46) !!!! call mxkrta(m,n) !!!! call xctmr1(46) !!!! call pipe_comparall(m,n, 'mxkrta, step') !!!! elseif (mxlktb) then !!!! call xctmr0(46) !!!! call mxkrtb(m,n) !!!! call xctmr1(46) !!!! call pipe_comparall(m,n, 'mxkrtb, step') !!!! else !!!! call xctmr0(46) !!!! call xctmr1(46) !!!! endif !mixed layer !!!! diagno = diagsv !!!! if (mxlkpp .or. mxlmy .or. mxlgiss) then !tsoff in mxkprf !!!! call xctmr0(47) !!!! call xctmr1(47) !!!! call xctmr0(48) !!!! call xctmr1(48) !!!! else ! mxlpwp has dypflg=2 !!!! call xctmr0(47) !!!! call convch(m,n) !!!! call xctmr1(47) !!!! call pipe_comparall(m,n, 'convch, step') !!!! if (dypflg.eq.1) then ! KPP-like, tsoff in diapf1 !!!! call xctmr0(48) !!!! call diapf1(m,n) !!!! call xctmr1(48) !!!! call pipe_comparall(m,n, 'diapf1, step') !!!! elseif (dypflg.eq.2) then ! explicit, tsoff in diapf2 !!!! call xctmr0(48) !!!! call diapf2(m,n) !!!! call xctmr1(48) !!!! call pipe_comparall(m,n, 'diapf2, step') !!!! else !!!! call xctmr0(48) !!!! call xctmr1(48) !!!! endif !!!! endif !diapycnal mixing !!!! if (incflg.ne.0 .and. incstp.gt.1) then !!!! call xctmr0(54) !!!! call incupd(n) !!!! call xctmr1(54) !!!! call pipe_comparall(m,n, 'incupd, step') !!!! else !!!! call xctmr0(54) !!!! call xctmr1(54) !!!! endif ! incremental update !!!! call xctmr0(49) !!!! call hybgen(m,n) !!!! call xctmr1(49) !!!! call pipe_comparall(m,n, 'hybgen, step') !!!! else ! isopyc !!!! call xctmr0(46) !!!! call mxkrtm(m,n) !!!! call xctmr1(46) !!!! call pipe_comparall(m,n, 'mxkrtm, step') !!!! call xctmr0(47) !!!! call convcm(m,n) !!!! call xctmr1(47) !!!! call pipe_comparall(m,n, 'convcm, step') !!!! call xctmr0(48) !!!! call diapf3(m,n) !!!! call xctmr1(48) !!!! call pipe_comparall(m,n, 'diapf3, step') !!!! call xctmr0(54) !!!! call xctmr1(54) !!!! call xctmr0(49) !!!! call xctmr1(49) !!!! endif !hybrid:isopyc !!!!c !!!!c --- update floats/moorings !!!! if (synflt) then !!!! nstepfl=nstepfl+1 !!!! if (nstepfl.eq.1) then !!!! iofl=1 !!!! call floats(m,n,0.0,iofl) !!!! endif !!!! if (mod(nstepfl,nfldta).eq.0) then !!!! iofl=0 !!!! if (mod(nstepfl,nflsam).eq.0) then !!!! iofl=1 !!!! endif !!!! call floats(m,n,real(rt%idd)+real(rt%ihh)/24.0,iofl) !!!! endif !!!! endif !synflt !!!!c !!!!#if defined (NERSC_VERSION) !!!!C --- Dump nesting conditions to inner models !!!! if (lnesto .and. imem==1) then !!!! call nestcondsave(dtime,n,nstep==nstep1+1) !!!! endif !!!!#endif /* NERSC */ !!!! !!!!c !!!!c --------------------------------------------------------------------------- !!!!c !!!!c --- output and diagnostic calculations !!!!c !!!!c --------------------------------------------------------------------------- !!!!c !!!!#if defined (NERSC_VERSION) !!!!c --- Update time info !!!! call year_day(dtime,refyear,rt,yrflag) !!!!#if defined (WEEKLY_AVERAGE) !!!!c --- For calculating/dumping average !!!! call ave_process(rt,m,n) !!!!#endif /* WEEKLY_AVERAGE */ !!!!c !!!!c --- Set flags loutput and outputday. Used to decide when to !!!!c --- dump restart/diagnostics, etc !!!! call set_outputday(rt,refyear,yrflag,baclin,nstep==nstep+1) !!!!c !!!!#if defined (DAILY_AVERAGE) !!!!c --- Process daily average (add/save if necessary) !!!! call daily_average(dtime,n,nstep==nstep1+1) !!!!#endif !!!!c --- Process grid points (add/save if necessary) !!!! call gridp_process(rt,dtime,m,n) !!!!c !!!!c --- Tide diag !!!! if (ltides.and.((mod(nstep,20)==0).or.(nstep==nstep+1))) then !!!! call tide_diag(n,dtime) !!!! if (mnproc==1) !!!! & write(lp,*) 'NERSC tides must be checked (new time input )' !!!! endif !!!!c !!!!c --- ice volume diag !!!!#if defined (ICE) !!!! call diag_icevol(dtime) !!!!#endif !!!!c !!!!c --- Status !!!! call prtstatus(dtime) !!!!c !!!!c --- If this is an input/output day ... !!!! if (loutput) then !!!! if (dday(outday)%res) then !!!! call save_restart_NERSC(nstep,dtime) !!!!Cc --- Save hycom and ice restart files !!!!C call save_restart_mem(rt,nstep,dtime) !!!!C call icemodels_save_restart(rt,nstep,dtime,nstep==nstep2) !!!!#if defined (WEEKLY_AVERAGE) !!!!c --- Save weekly average restart file !!!! call save_restart_ave(rt) !!!!c !!!!#endif /* WEEKLY_AVERAGE */ !!!!Cc --- Save random forcing !!!!C if (randf) call save_restart_rand_ab(rt) !!!!Cc !!!!C#if defined (SINGLE_RESTART) /* We dump ssh when this flag is set - its convenient*/ !!!!Cc --- Also put ssh field in restart file - useful for assimilation !!!!C call restart_appendfield(trim(restartfile(rt)), !!!!C & 'ssh ',0,1,srfhgt ,ierr) !!!!Cc !!!!C#endif !!!! end if !!!! end if !!!!#endif /* NERSC_VERSION */ !!!!c !!!! lfatal = .false. !!!! diag_tide = tidflg.gt.0 .and. !!!! & mod(dtime+dsmall,hours1).lt.dsmall2 !at least hourly !!!! if (diagno .or. histry .or. hitile .or. hisurf .or. histmn .or. !!!! & nstep.eq.nstep0+1 .or. !!!! & diag_tide .or. !!!! & mod(dtime+dsmall,ddsurf).lt.dsmall2 .or. !!!! & mod(dtime+dsmall, days1).lt.dsmall2 ) then ! at least daily !!!!c !!!!#if defined(USE_CCSM3) !!!! jday = days_in_prior_months(imonth)+iday !!!!#else !!!! call forday(dtime,yrflag, iyear,jday,ihour) !!!!#endif !!!! write(c_ydh,'('' ('',i4.4,''/'',i3.3,1x,i2.2,'')'')') !!!! & iyear,jday,ihour !!!!c !!!!c --- diagnose mean sea surface height !!!! !!!! if (.not.allocated(sminy)) then !!!! allocate( sminy(1:jj), smaxy(1:jj) ) !!!! endif !!!!!$OMP PARALLEL DO PRIVATE(j,l,i) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1,jj !!!! sminy(j)= huge !!!! smaxy(j)=-huge !!!! do l=1,isp(j) !!!! do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) !!!! if (tidflg.gt.0) then !!!! util2(i,j)=(srfhgt(i,j)/g)**2*scp2(i,j) !!!! endif !!!! util3(i,j)=srfhgt(i,j)*scp2(i,j) !!!! util4(i,j)=montg1(i,j)*scp2(i,j) !!!! sminy(j)=min(sminy(j),srfhgt(i,j)) !!!! smaxy(j)=max(smaxy(j),srfhgt(i,j)) !!!! enddo !!!! enddo !!!! enddo !!!!!$OMP END PARALLEL DO !!!! smin=minval(sminy(1:jj)) !!!! smax=maxval(smaxy(1:jj)) !!!! call xcminr(smin) !!!! call xcmaxr(smax) !!!! call xcsum( dsum, util3,ip) !!!! call xcsum( dsmt, util4,ip) !!!! sum=dsum !!!! smt=dsmt !!!! if (mnproc.eq.1) then !!!! write (lp,'(i9,a, !!!! & '' mean SSH (mm):'',f8.2, !!!! & '' ('',1pe8.1,'' to '',e8.1,'')'')') !!!! & nstep,c_ydh, !!!! & sum/(area*thref*onemm),smin/(thref*onemm),smax/(thref*onemm) !!!!* write (lp,'(i9,a, !!!!* . '' mean MontgPot (mm):'',f8.2)') !!!!* . nstep,c_ydh, !!!!* . smt/(area*thref*onemm) !!!! call flush(lp) !!!! write(nod,'(i9,a, !!!! & '' mean SSH (mm):'',f8.2, !!!! & '' ('',1pe8.1,'' to '',e8.1,'')'')') !!!! & nstep,c_ydh, !!!! & sum/(area*thref*onemm),smin/(thref*onemm),smax/(thref*onemm) !!!!* write(nod,'(i9,a, !!!!* . '' mean MontgPot (mm):'',f8.2)') !!!!* . nstep,c_ydh, !!!!* . smt/(area*thref*onemm) !!!! call flush(nod) !!!! endif !1st tile !!!!c --- NaN detection. !!!! if (hycom_isnaninf(smin) .or. !!!! & hycom_isnaninf(sum) .or. !!!! & hycom_isnaninf(smax) .or. !!!! & hycom_isnaninf(smt) ) then !!!! if (mnproc.eq.1) then !!!! write(lp,*) !!!! write(lp,*) 'error - NaN or Inf detected' !!!! write(lp,*) !!!! call flush(lp) !!!! endif !1st tile !!!! lfatal = .true. !delay exit to allow archive output !!!! endif !NaN !!!!c !!!! if (tidflg.gt.0) then !!!! call xcsum( dsms, util2,ip) !!!! sms=dsms/area !!!!c !!!! call xctilr(u( 1-nbdy,1-nbdy,1,n),1,kk, 1,1, halo_uv) !!!! call xctilr(v( 1-nbdy,1-nbdy,1,n),1,kk, 1,1, halo_vv) !!!! call xctilr(ubavg(1-nbdy,1-nbdy, n),1, 1, 1,1, halo_uv) !!!! call xctilr(vbavg(1-nbdy,1-nbdy, n),1, 1, 1,1, halo_vv) !!!!c !!!! dskea=0.0d0 !!!! do k=1,kk !!!!!$OMP PARALLEL DO PRIVATE(j,l,i,utotp,vtotp) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1,jj !!!! do l=1,isp(j) !!!! do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) !!!! utotp=0.5*( u(i, j,k,n)+ubavg(i, j,n) + !!!! & u(i+1,j,k,n)+ubavg(i+1,j,n) ) !!!! vtotp=0.5*( v(i,j, k,n)+vbavg(i,j, n) + !!!! & v(i,j+1,k,n)+vbavg(i,j+1,n) ) !!!! util4(i,j)=dp(i,j,k,n)*scp2(i,j)* !!!! & 0.5*(1000.0+th3d(i,j,k,n)+thbase)* !!!! & (utotp**2+vtotp**2) !!!! enddo !!!! enddo !!!! enddo !!!!!$OMP END PARALLEL DO !!!! call xcsum(dske, util4,ip) !!!! dskea=dskea+dske !!!! enddo !k !!!! sum=dskea/(area*onem) !!!! if (mnproc.eq.1) then !!!! write (lp,'(i9,a, !!!! & '' region-wide mean KE: '',f20.10)') !!!! & nstep,c_ydh, !!!! & sum !!!! write (lp,'(i9,a, !!!! & '' region-wide mean APE:'',f20.10)') !!!! & nstep,c_ydh, !!!! & sms*0.5*g*(1000.0+thbase) !!!! endif !!!! endif !tidflg !!!!* else !!!!* if (mnproc.eq.1) then !!!!* write (lp,'('' time step ='',i9)') nstep !!!!* call flush(lp) !!!!* endif !1st tile !!!! endif !daily or hourly !!!!c !!!!c --- diagnose heat/salt flux, ice, layer thickness and temperature, !!!!c --- mean temperature and mean kinetic energy !!!!c --- note that mixed-layer fields must be switched on in mxkprf !!!! if (diagno .or. !!!! & nstep.eq.nstep0+1 .or. !!!! & mod(dtime+dsmall,days6).lt.dsmall2) then ! at least every 6 days !!!!c !!!!#if defined(USE_CCSM3) !!!! jday = days_in_prior_months(imonth)+iday !!!!#else !!!! call forday(dtime,yrflag, iyear,jday,ihour) !!!!#endif !!!! write(c_ydh,'('' ('',i4.4,''/'',i3.3,1x,i2.2,'')'')') !!!! & iyear,jday,ihour !!!!c !!!! if (thermo .or. sstflg.gt.0 .or. srelax) then !!!!!$OMP PARALLEL DO PRIVATE(j,l,i) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1,jj !!!! do l=1,isp(j) !!!! do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) !!!! util1(i,j)=buoflx(i,j)*scp2(i,j) !!!! util2(i,j)=bhtflx(i,j)*scp2(i,j) !!!! enddo !!!! enddo !!!! enddo !!!!!$OMP END PARALLEL DO !!!! call xcsum(dsum, util1,ip) !!!! call xcsum(dsmt, util2,ip) !!!! sum= (dsum*1.00D9)/area ! 1.e9*m**2/sec**3 !!!! smt= (dsmt*1.00D9)/area ! 1.e9*m**2/sec**3 !!!! if (mnproc.eq.1) then !!!! write (lp, '(i9,a, !!!! & '' mean BFL (m^2/s^3):'',f8.2, !!!! & '' hfl:'',f8.2)') !!!! & nstep,c_ydh, !!!! & sum,smt !!!! call flush(lp) !!!! write (nod,'(i9,a, !!!! & '' mean BFL (m^2/s^3):'',f8.2, !!!! & '' hfl:'',f8.2)') !!!! & nstep,c_ydh, !!!! & sum,smt !!!! call flush(nod) !!!! endif !1st tile !!!!c !!!!!$OMP PARALLEL DO PRIVATE(j,l,i) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1,jj !!!! do l=1,isp(j) !!!! do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) !!!! util1(i,j)=surflx(i,j)*scp2(i,j) !!!! util2(i,j)=mixflx(i,j)*scp2(i,j) !!!! util3(i,j)=sstflx(i,j)*scp2(i,j) !!!! enddo !!!! enddo !!!! enddo !!!!!$OMP END PARALLEL DO !!!! call xcsum(dsum, util1,ip) !!!! call xcsum(dsmt, util2,ip) !!!! call xcsum(d3, util3,ip) !!!! sum= dsum/area !!!! smt= dsmt/area !!!! smr= d3 /area !!!! if (mnproc.eq.1) then !!!! write (lp, '(i9,a, !!!! & '' mean HFLUX (w/m^2):'',f8.2, !!!! & '' sst:'',f8.2, !!!! & '' ml:'',f8.2)') !!!! & nstep,c_ydh, !!!! & sum,smr,smt !!!! call flush(lp) !!!! write (nod,'(i9,a, !!!! & '' mean HFLUX (w/m^2):'',f8.2, !!!! & '' sst:'',f8.2, !!!! & '' ml:'',f8.2)') !!!! & nstep,c_ydh, !!!! & sum,smr,smt !!!! call flush(nod) !!!! endif !1st tile !!!!c !!!!!$OMP PARALLEL DO PRIVATE(j,l,i) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1,jj !!!! do l=1,isp(j) !!!! do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) !!!! util1(i,j)=salflx(i,j)*scp2(i,j)/saln(i,j,1,n) !!!! util2(i,j)=sssflx(i,j)*scp2(i,j)/saln(i,j,1,n) !!!! enddo !!!! enddo !!!! enddo !!!!!$OMP END PARALLEL DO !!!! call xcsum(dsms, util1,ip) !!!! call xcsum(dsum, util2,ip) !!!! sms=-(dsms*thref*7.0D0*8.64D7)/area ! P-E in mm/week !!!! smr=-(dsum*thref*7.0D0*8.64D7)/area ! P-E in mm/week !!!! if (mnproc.eq.1) then !!!! write (lp, '(i9,a, !!!! & '' mean WFLUX (mm/wk):'',f8.2, !!!! & '' sss:'',f8.2)') !!!! & nstep,c_ydh, !!!! & sms,smr !!!! call flush(lp) !!!! write (nod,'(i9,a, !!!! & '' mean WFLUX (mm/wk):'',f8.2, !!!! & '' sss:'',f8.2)') !!!! & nstep,c_ydh, !!!! & sms,smr !!!! call flush(nod) !!!! endif !1st tile !!!! endif !!!!c !!!! if (icegln) then ! basin-wide ice !!!!!$OMP PARALLEL DO PRIVATE(j,l,i) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1,jj !!!! do l=1,isp(j) !!!! do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) !!!! if (covice(i,j).ne.0.0) then !!!! util2(i,j)=covice(i,j)* scp2(i,j) !!!! util3(i,j)= thkice(i,j)*scp2(i,j) !!!! util4(i,j)=covice(i,j)*temice(i,j)*scp2(i,j) !!!! else !!!! util2(i,j)=0.0 !!!! util3(i,j)=0.0 !!!! util4(i,j)=0.0 !!!! endif !!!! enddo !!!! enddo !!!! enddo !!!!!$OMP END PARALLEL DO !!!! call xcsum(d2, util2,ip) !!!! call xcsum(d3, util3,ip) !!!! call xcsum(d4, util4,ip) !!!! if (d2.gt.0.0d0) then !!!! sum=d3/d2 !average ice thickness, where there is ice !!!! smt=d4/d2 !average ice temperature, where there is ice !!!! sms=d2/area * 100.0 !ice coverage, percent of total area !!!! else !!!! sum=0.0 !!!! smt=0.0 !!!! sms=0.0 !!!! endif !!!! if (mnproc.eq.1) then !!!! write (lp,'(i9,a, !!!! & '' mean ice thk. (m):'',f8.2, !!!! & '' temp:'',f7.3, !!!! & '' pcen:'',f7.3)') !!!! & nstep,c_ydh, !!!! & sum,smt,sms !!!! call flush(lp) !!!! write(nod,'(i9,a, !!!! & '' mean ice thk. (m):'',f8.2, !!!! & '' temp:'',f7.3, !!!! & '' pcen:'',f7.3)') !!!! & nstep,c_ydh, !!!! & sum,smt,sms !!!! call flush(nod) !!!! endif !1st tile !!!! endif ! icegln !!!!c !!!! if (nreg.ne.0 .and. icegln) then ! southern hemisphere ice !!!! d2a = d2 !!!! d3a = d3 !!!! d4a = d4 !!!!!$OMP PARALLEL DO PRIVATE(j,l,i) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1,jj !!!! do l=1,isp(j) !!!! do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) !!!! if (covice(i,j).ne.0.0 .and. !!!! & plat(i,j).lt.0.0 ) then !!!! util2(i,j)=covice(i,j)* scp2(i,j) !!!! util3(i,j)= thkice(i,j)*scp2(i,j) !!!! util4(i,j)=covice(i,j)*temice(i,j)*scp2(i,j) !!!! else !!!! util2(i,j)=0.0 !!!! util3(i,j)=0.0 !!!! util4(i,j)=0.0 !!!! endif !!!! enddo !!!! enddo !!!! enddo !!!!!$OMP END PARALLEL DO !!!! call xcsum(d2, util2,ip) !!!! call xcsum(d3, util3,ip) !!!! call xcsum(d4, util4,ip) !!!! if (d2.gt.0.0d0) then !!!! sum=d3/d2 !average ice thickness, where there is ice in S.H. !!!! smt=d4/d2 !average ice temperature, where there is ice in S.H. !!!! sms=d2/area * 100.0 !S.H. ice coverage, percent of total area !!!! else !!!! sum=0.0 !!!! smt=0.0 !!!! sms=0.0 !!!! endif !!!! if (mnproc.eq.1) then !!!! write (lp,'(i9,a, !!!! & '' mean SH I thk. (m):'',f8.2, !!!! & '' temp:'',f7.3, !!!! & '' pcen:'',f7.3)') !!!! & nstep,c_ydh, !!!! & sum,smt,sms !!!! call flush(lp) !!!! write(nod,'(i9,a, !!!! & '' mean SH I thk. (m):'',f8.2, !!!! & '' temp:'',f7.3, !!!! & '' pcen:'',f7.3)') !!!! & nstep,c_ydh, !!!! & sum,smt,sms !!!! call flush(nod) !!!! endif !1st tile !!!!c !!!! d2 = d2a - d2 !!!! d3 = d3a - d3 !!!! d4 = d4a - d4 !!!! if (d2.gt.0.0d0) then !!!! sum=d3/d2 !average ice thickness, where there is ice in N.H. !!!! smt=d4/d2 !average ice temperature, where there is ice in N.H. !!!! sms=d2/area * 100.0 !N.H. ice coverage, percent of total area !!!! else !!!! sum=0.0 !!!! smt=0.0 !!!! sms=0.0 !!!! endif !!!! if (mnproc.eq.1) then !!!! write (lp,'(i9,a, !!!! & '' mean NH I thk. (m):'',f8.2, !!!! & '' temp:'',f7.3, !!!! & '' pcen:'',f7.3)') !!!! & nstep,c_ydh, !!!! & sum,smt,sms !!!! call flush(lp) !!!! write(nod,'(i9,a, !!!! & '' mean NH I thk. (m):'',f8.2, !!!! & '' temp:'',f7.3, !!!! & '' pcen:'',f7.3)') !!!! & nstep,c_ydh, !!!! & sum,smt,sms !!!! call flush(nod) !!!! endif !1st tile !!!! endif ! icegln .and. nreg.ne.0 !!!!c !!!! if (icegln .and. !!!! & (icmflg.eq.2 .or. ticegr.eq.0.0 .or. !!!! & lwflag.eq.2 .or. sstflg.eq.2 )) then !!!!c !!!!c --- ice mask !!!!c !!!!!$OMP PARALLEL DO PRIVATE(j,l,i,tsur) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1,jj !!!! do l=1,isp(j) !!!! do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) !!!!c !!!!c --- don't allow a new tsur maximum, to preserve sea ice !!!!c !!!! if (yrflag.lt.2) then !!!! tsur = min( max( surtmp(i,j,l0), surtmp(i,j,l1), !!!! & surtmp(i,j,l2), surtmp(i,j,l3) ), !!!! & surtmp(i,j,l0)*w0+surtmp(i,j,l1)*w1+ !!!! & surtmp(i,j,l2)*w2+surtmp(i,j,l3)*w3 ) !!!! else !!!! tsur = min( max( surtmp(i,j,l0), surtmp(i,j,l1) ), !!!! & surtmp(i,j,l0)*w0+surtmp(i,j,l1)*w1 ) !!!! endif !!!! if (tsur.le.tfrz_n) then !!!! util1(i,j) = scp2(i,j) !!!! if (plat(i,j).lt.0.0) then !!!! util2(i,j) = scp2(i,j) !!!! else !!!! util2(i,j) = 0.0 !!!! endif !!!! else !!!! util1(i,j) = 0.0 !!!! util2(i,j) = 0.0 !!!! endif !!!! enddo !!!! enddo !!!! enddo !!!!!$OMP END PARALLEL DO !!!! call xcsum(dsmt, util1,ip) !!!! call xcsum(dsms, util2,ip) !!!! smt=dsmt/area * 100.0 !!!! sms=dsms/area * 100.0 !!!! if (mnproc.eq.1) then !!!! write (lp,'(i9,a, !!!! & '' mean ice mask pcen:'',f9.3, !!!! & '' S.H:'',f7.3, !!!! & '' N.H:'',f7.3)') !!!! & nstep,c_ydh, !!!! & smt,sms,smt-sms !!!! call flush(lp) !!!! write(nod,'(i9,a, !!!! & '' mean ice mask pcen:'',f9.3, !!!! & '' S.H:'',f7.3, !!!! & '' N.H:'',f7.3)') !!!! & nstep,c_ydh, !!!! & smt,sms,smt-sms !!!! call flush(nod) !!!! endif !1st tile !!!! endif !ice mask statistics !!!!c !!!!!$OMP PARALLEL DO PRIVATE(j,l,i) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1,jj !!!! do l=1,isp(j) !!!! do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) !!!! util1(i,j)=dpmixl(i,j,n)*scp2(i,j) !!!! util2(i,j)=dpmixl(i,j,n)*scp2(i,j)*tmix(i,j) !!!! util3(i,j)=dpmixl(i,j,n)*scp2(i,j)*smix(i,j) !!!! enddo !!!! enddo !!!! enddo !!!!!$OMP END PARALLEL DO !!!! call xcsum(dsum, util1,ip) !!!! call xcsum(dsmt, util2,ip) !!!! call xcsum(dsms, util3,ip) !!!! if (dsum.ne.0.0d0) then !!!! sum=dsum/(area*onem) !!!! smt=dsmt/dsum !!!! sms=dsms/dsum !!!! else !!!! sum=0.0 !!!! smt=0.0 !!!! sms=0.0 !!!! endif !!!! if (mnproc.eq.1) then !!!! write (lp,'(i9,a, !!!! & '' mean mixl thk. (m):'',f8.2, !!!! & '' temp:'',f7.3, !!!! & '' saln:'',f7.3)') !!!! & nstep,c_ydh, !!!! & sum,smt,sms !!!! call flush(lp) !!!! write(nod,'(i9,a, !!!! & '' mean mixl thk. (m):'',f8.2, !!!! & '' temp:'',f7.3, !!!! & '' saln:'',f7.3)') !!!! & nstep,c_ydh, !!!! & sum,smt,sms !!!! call flush(nod) !!!! endif !1st tile !!!!c !!!! if (relaxf .or. sstflg.ne.0) then !!!!c !!!!c --- mean surface climatology. !!!!c !!!!!$OMP PARALLEL DO PRIVATE(j,l,i) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1,jj !!!! do l=1,isp(j) !!!! do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) !!!! if (relaxf .and. sstflg.le.1) then !!!! util1(i,j)=scp2(i,j)* !!!! & (twall(i,j,1,lc0)*wc0+twall(i,j,1,lc1)*wc1 !!!! & +twall(i,j,1,lc2)*wc2+twall(i,j,1,lc3)*wc3) !!!! else !synoptic observed sst !!!! util1(i,j)=scp2(i,j)* !!!! & (seatmp(i,j,l0)*w0+seatmp(i,j,l1)*w1 !!!! & +seatmp(i,j,l2)*w2+seatmp(i,j,l3)*w3) !!!! endif !!!! if (relaxf) then !sss !!!! util2(i,j)=scp2(i,j)* !!!! & (swall(i,j,1,lc0)*wc0+swall(i,j,1,lc1)*wc1 !!!! & +swall(i,j,1,lc2)*wc2+swall(i,j,1,lc3)*wc3) !!!! else !synoptic observed surface temperature !!!! util1(i,j)=scp2(i,j)* !!!! & (surtmp(i,j,l0)*w0+surtmp(i,j,l1)*w1 !!!! & +surtmp(i,j,l2)*w2+surtmp(i,j,l3)*w3) !!!! endif !!!! enddo !!!! enddo !!!! enddo !!!!!$OMP END PARALLEL DO !!!! call xcsum(dsmt, util1,ip) !!!! call xcsum(dsms, util2,ip) !!!! smt=dsmt/area !!!! sms=dsms/area !!!! if (mnproc.eq.1) then !!!! if (relaxf) then !!!! write (lp,'(i9,a, !!!! & '' mean clim thk. (m):'',f8.2, !!!! & '' sst:'',f7.3, !!!! & '' sss:'',f7.3)') !!!! & nstep,c_ydh, !!!! & thkmin,smt,sms !!!! call flush(lp) !!!! write(nod,'(i9,a, !!!! & '' mean clim thk. (m):'',f8.2, !!!! & '' sst:'',f7.3, !!!! & '' sss:'',f7.3)') !!!! & nstep,c_ydh, !!!! & thkmin,smt,sms !!!! call flush(nod) !!!! else !.not.relaxf !!!! write (lp,'(i9,a, !!!! & '' mean clim thk. (m):'',f8.2, !!!! & '' sst:'',f7.3, !!!! & '' surt:'',f7.3)') !!!! & nstep,c_ydh, !!!! & thkmin,smt,sms !!!! call flush(lp) !!!! write(nod,'(i9,a, !!!! & '' mean clim thk. (m):'',f8.2, !!!! & '' sst:'',f7.3, !!!! & '' surt:'',f7.3)') !!!! & nstep,c_ydh, !!!! & thkmin,smt,sms !!!! call flush(nod) !!!! endif !relaxf:else !!!! endif !1st tile !!!! endif !!!!c !!!! call xctilr(u( 1-nbdy,1-nbdy,1,n),1,kk, 1,1, halo_uv) !!!! call xctilr(v( 1-nbdy,1-nbdy,1,n),1,kk, 1,1, halo_vv) !!!! call xctilr(ubavg(1-nbdy,1-nbdy, n),1, 1, 1,1, halo_uv) !!!! call xctilr(vbavg(1-nbdy,1-nbdy, n),1, 1, 1,1, halo_vv) !!!!c !!!! dsuma=0.0d0 !!!! dsmta=0.0d0 !!!! dsmsa=0.0d0 !!!! dsmra=0.0d0 !!!! dskea=0.0d0 !!!! do k=1,kk !!!!!$OMP PARALLEL DO PRIVATE(j,l,i,utotp,vtotp) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1,jj !!!! do l=1,isp(j) !!!! do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) !!!! utotp=0.5*( u(i, j,k,n)+ubavg(i, j,n) + !!!! & u(i+1,j,k,n)+ubavg(i+1,j,n) ) !!!! vtotp=0.5*( v(i,j, k,n)+vbavg(i,j, n) + !!!! & v(i,j+1,k,n)+vbavg(i,j+1,n) ) !!!! util1(i,j)=dp(i,j,k,n)*scp2(i,j) !!!! util2(i,j)=dp(i,j,k,n)*scp2(i,j)*temp(i,j,k,n) !!!! util3(i,j)=dp(i,j,k,n)*scp2(i,j)*saln(i,j,k,n) !!!! util5(i,j)=dp(i,j,k,n)*scp2(i,j)*th3d(i,j,k,n) !!!! util4(i,j)=dp(i,j,k,n)*scp2(i,j)* !!!! & 0.5*(1000.0+th3d(i,j,k,n)+thbase)* !!!! & (utotp**2+vtotp**2) !!!! enddo !!!! enddo !!!! enddo !!!!!$OMP END PARALLEL DO !!!! call xcsum(dsum, util1,ip) !!!! call xcsum(dsmt, util2,ip) !!!! call xcsum(dsms, util3,ip) !!!! call xcsum(dsmr, util5,ip) !!!! call xcsum(dske, util4,ip) !!!! dsuma=dsuma+dsum !!!! dsmta=dsmta+dsmt !!!! dsmsa=dsmsa+dsms !!!! dsmra=dsmra+dsmr !!!! dskea=dskea+dske !!!! if (dsum.ne.0.0d0) then !!!! sum=dsum/(area*onem) !!!! smt=dsmt/dsum !!!! sms=dsms/dsum !!!! else !!!! sum=0.0 !!!! smt=0.0 !!!! sms=0.0 !!!! endif !!!! if (mnproc.eq.1) then !!!! write (lp,'(i9,a, !!!! & '' mean L '',i2,'' thk. (m):'',f8.2, !!!! & '' temp:'',f7.3, !!!! & '' saln:'',f7.3)') !!!! & nstep,c_ydh, !!!! & k,sum,smt,sms !!!! call flush(lp) !!!! write(nod,'(i9,a, !!!! & '' mean L '',i2,'' thk. (m):'',f8.2, !!!! & '' temp:'',f7.3, !!!! & '' saln:'',f7.3)') !!!! & nstep,c_ydh, !!!! & k,sum,smt,sms !!!! call flush(nod) !!!! endif !1st tile !!!! enddo !k !!!! sum=dskea/(area*onem) !!!! smt=dsmta/dsuma !!!! sms=dsmsa/dsuma !!!! smr=dsmra/dsuma !!!! if (mnproc.eq.1) then !!!! write (lp,'(i9,a, !!!! & '' region-wide mean Kin. Energy:'',f20.10)') !!!! & nstep,c_ydh, !!!! & sum !!!! write (lp,'(i9,a, !!!! & '' region-wide mean Temperature:'',f20.10)') !!!! & nstep,c_ydh, !!!! & smt !!!! write (lp,'(i9,a, !!!! & '' region-wide mean Salinity: '',f20.10)') !!!! & nstep,c_ydh, !!!! & sms !!!! write (lp,'(i9,a, !!!! & '' region-wide mean Density Dev:'',f20.10)') !!!! & nstep,c_ydh, !!!! & smr !!!! call flush(lp) !!!! write(nod,'(i9,a, !!!! & '' region-wide mean Kin. Energy:'',f20.10)') !!!! & nstep,c_ydh, !!!! & sum !!!! write(nod,'(i9,a, !!!! & '' region-wide mean Temperature:'',f20.10)') !!!! & nstep,c_ydh, !!!! & smt !!!! write(nod,'(i9,a, !!!! & '' region-wide mean Salinity: '',f20.10)') !!!! & nstep,c_ydh, !!!! & sms !!!! write(nod,'(i9,a, !!!! & '' region-wide mean Density Dev:'',f20.10)') !!!! & nstep,c_ydh, !!!! & smr !!!! call flush(nod) !!!! endif !1st tile !!!!c !!!! do ktr= 1,ntracr !!!! dsumtr(ktr)=0.0d0 !!!! do k=1,kk !!!!!$OMP PARALLEL DO PRIVATE(j,l,i) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1,jj !!!! do l=1,isp(j) !!!! do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) !!!! util2(i,j)=dp(i,j,k,n)*scp2(i,j)*tracer(i,j,k,n,ktr) !!!! enddo !i !!!! enddo !l !!!! enddo !j !!!!!$OMP END PARALLEL DO !!!! call xcsum(dsmt, util2,ip) !!!! dsumtr(ktr)=dsumtr(ktr)+dsmt !!!! enddo !k !!!! smt=dsumtr(ktr)/dsuma !dsuma still good from K.E. loops !!!! if (mnproc.eq.1) then !!!! write (lp,'(i9,a, !!!! & '' region-wide mean tracer'',i3.2, !!!! & '': '',f20.10)') !!!! & nstep,c_ydh, !!!! & ktr,smt !!!! call flush(lp) !!!! write(nod,'(i9,a, !!!! & '' region-wide mean tracer'',i3.2, !!!! & '': '',f20.10)') !!!! & nstep,c_ydh, !!!! & ktr,smt !!!! call flush(nod) !!!! endif !1st tile !!!!c --- NaN detection, for each tracer. !!!! if (hycom_isnaninf(smt)) then !!!! if (mnproc.eq.1) then !!!! write(lp,*) !!!! write(lp,*) 'error - NaN or Inf detected' !!!! write(lp,*) !!!! call flush(lp) !!!! endif !1st tile !!!! lfatal = .true. !delay exit to allow archive output !!!! endif !NaN !!!! if (ktr.ge.3 .and. trcflg(ktr-2).eq.903) then !NPZ !!!! smt=(dsumtr(ktr)+dsumtr(ktr-1)+dsumtr(ktr-2))/dsuma !!!! if (mnproc.eq.1) then !!!! write (lp,'(i9,a, !!!! & '' region-wide mean N+P+Z: '',f20.10)') !!!! & nstep,c_ydh, !!!! & smt !!!! call flush(lp) !!!! write(nod,'(i9,a, !!!! & '' region-wide mean N+P+Z: '',f20.10)') !!!! & nstep,c_ydh, !!!! & smt !!!! call flush(nod) !!!! endif !1st tile !!!! elseif (ktr.ge.4 .and. trcflg(ktr-3).eq.904) then !NPZD !!!! smt=(dsumtr(ktr) +dsumtr(ktr-1)+ !!!! & dsumtr(ktr-2)+dsumtr(ktr-3))/dsuma !!!! if (mnproc.eq.1) then !!!! write (lp,'(i9,a, !!!! & '' region-wide mean N+P+Z+D: '',f20.10)') !!!! & nstep,c_ydh, !!!! & smt !!!! call flush(lp) !!!! write(nod,'(i9,a, !!!! & '' region-wide mean N+P+Z+D: '',f20.10)') !!!! & nstep,c_ydh, !!!! & smt !!!! call flush(nod) !!!! endif !1st tile !!!! endif !NPZ:NPZD !!!! enddo !ktr !!!! endif !diagno ... !!!!c !!!!c --- diagnose meridional overturning and heat flux !!!! if (mod(dtime+dsmall,dmonth).lt.dsmall2 .and. imem==1) then !!!! call xctmr0(52) !!!! call overtn(dtime,dyear) !!!! call xctmr1(52) !!!! elseif (nstep.ge.nstep2) then !!!! call xctmr0(52) !!!! call overtn(dtime,dyear) !!!! call xctmr1(52) !!!! endif !!!!c !!!! if (meanfq.ne.0.0) then !!!! if (.not. histmn) then !!!! if (tidflg.gt.0 .and. meanfq.eq.1.0) then !!!! if (mod(dtime+dsmall,hours1).lt.dsmall2) then !!!! call mean_add(n, 1.0) !25-hour average daily !!!! endif !!!! else !!!! call mean_add(n, 1.0) !!!! endif !!!! else ! histmn !!!! if (tidflg.gt.0 .and. meanfq.eq.1.0) then !!!! call mean_add(n, 1.0) !25-hour average daily !!!! else !!!! call mean_add(n, 0.5) !!!! endif !!!!c !!!! call xctmr0(53) !!!!c !!!!c --- output to mean archive file !!!!c !!!! call mean_end(dtime) !!!!#if defined(USE_CCSM3) !!!! jday = days_in_prior_months(imonth)+iday !!!!#else !!!! call forday(time_ave,yrflag, iyear,jday,ihour) !!!!#endif !!!! call mean_archiv(n, iyear,jday,ihour) !!!!c !!!! call mean_zero(dtime) !!!! call mean_add(n, 0.5) !!!!c !!!! call xctmr1(53) !!!! endif ! histmn !!!! endif !meanfq !!!!c !!!! if (histry .or. hitile .or. hisurf .or. lfatal) then !!!! call xctmr0(53) !!!!c !!!!c --- output to archive file !!!!c !!!!#if defined(USE_CCSM3) !!!! jday = days_in_prior_months(imonth)+iday !!!!#else !!!! call forday(dtime,yrflag, iyear,jday,ihour) !!!!#endif !!!!c !!!!!$OMP PARALLEL DO PRIVATE(j,k,l,i) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1,jj !!!! if (isopyc .or. mxlkrt) then !!!! do l=1,isu(j) !!!! do i=max(1,ifu(j,l)),min(ii,ilu(j,l)) !!!! umix(i,j)=u(i,j,1,n) !!!! enddo !!!! enddo !!!! do l=1,isv(j) !!!! do i=max(1,ifv(j,l)),min(ii,ilv(j,l)) !!!! vmix(i,j)=v(i,j,1,n) !!!! enddo !!!! enddo !!!! endif !!!! if (histry) then !!!! do k= 1,kk !!!! do l=1,isp(j) !!!! do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) !!!!c --- convert diapycnal thickness changes into !!!!c --- actual interface fluxes !!!! if (k.gt.1) then !!!! diaflx(i,j,k)=diaflx(i,j,k)/(2.*onem) + !!!! & diaflx(i,j,k-1) !!!! else !!!! diaflx(i,j,k)=diaflx(i,j,k)/(2.*onem) !!!! endif !!!! enddo !!!! enddo !!!! enddo !k !!!! endif !!!! enddo !j !!!!!$OMP END PARALLEL DO !!!!c !!!! if (lfatal) then !write archive and exit !!!! if (mnproc.eq.1) then !!!! write (intvl,'(i3.3)') 0 !!!! endif !1st tile !!!! call archiv(n, kk, iyear,jday,ihour, intvl) !!!! call xcstop('(hycom)') !!!! stop '(hycom)' !won't get here !!!! else !!!! if (mnproc.eq.1) then !!!! write (intvl,'(i3.3)') int(dtime-timav+dsmall) !!!! endif !1st tile !!!! if (hisurf .and. .not. histry) then !!!! call archiv(n, 1, iyear,jday,ihour, intvl) !!!! elseif (histry) then !!!! call archiv(n, kk, iyear,jday,ihour, intvl) !!!! endif !hisurf:histry !!!! if (hitile) then !!!! call archiv_tile(n, kk, iyear,jday,ihour, intvl) !!!! endif !!!! endif !lfatal:else !!!!c !!!! if (histry) then !!!!!$OMP PARALLEL DO PRIVATE(j,k,l,i) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1,jj !!!! do k= 1,kk !!!! do l=1,isp(j) !!!! do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) !!!! diaflx(i,j,k)=0.0 !!!! enddo !!!! enddo !!!! enddo !!!! enddo !!!!!$OMP END PARALLEL DO !!!!c !!!! timav=time !!!! endif !!!! call xctmr1(53) !!!! endif ! histry.or.hitile.or.hisurf.or.lfatal !!!!c !!!!#if defined(USE_ESMF) !!!! if (put_export) then !!!!c !!!!c --- fill the ESMF Export State. !!!!c !!!! call Export_ESMF !!!! call forday(dtime,yrflag, iyear,jday,ihour) !!!! call Archive_ESMF(iyear,jday,ihour) !!!! endif !!!!#endif !!!!c !!!! if (restrt) then !!!! call xctmr0(51) !!!!c !!!!c --- output to restart and flux statitics files !!!!c !!!! if (mxlkpp) then !!!!!$OMP PARALLEL DO PRIVATE(j,l,i) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1,jj !!!! do l=1,isp(j) !!!! do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) !!!! dpmixl(i,j,m) = dpmixl(i,j,n) !!!! enddo !!!! enddo !!!! enddo !!!!!$OMP END PARALLEL DO !!!! endif !!!!c !!!!#if defined(USE_CCSM3) !!!! call ccsm3_time_date_stamp (ccsm3_string, 'ymds') !!!! !!!! flnmra = trim(flnmrsod)// '.'//trim(ccsm3_string) !!!! flnmrb = trim(flnmrsod)//'h.'//trim(ccsm3_string) !!!! !!!! if (mnproc == 1) then !!!! open(1,file=trim(pointer_filename),form='formatted', !!!! & status='unknown') !!!! write(1,'(a)')trim(flnmra) !!!! write(1,'(a)')trim(flnmrb) !!!! write(1,*)time ! real, in days !!!! write(1,*)dtime ! real8, in days !!!! write(1,*)nstep,iyear,imonth,iday ! integer !!!! &, elapsed_days !!!! write(1,*)eom, eoy ! logical !!!! close(1) !!!!c !!!! write(lp,'(2a)')'RESTART: (mainloop) to file: ',trim(flnmra) !!!! write(lp,*) ' time = ',time !!!! write(lp,*) ' nstep,iyear,imonth,iday : ', !!!! & nstep,iyear,imonth,iday !!!! call flush(lp) !!!! endif !mnproc == 1 !!!!#else !!!! call forday(dtime,yrflag, iyear,jday,ihour) !!!! if (mnproc.eq.1) then !!!! write (lp,'(a,i9, 9x,a,i6.4, 9x,a,i5.3, 9x,a,i4.2)') !!!! & ' time step',nstep, !!!! & 'y e a r', iyear, !!!! & 'd a y', jday, !!!! & 'h o u r', ihour !!!! call flush(lp) !!!! endif !1st tile !!!!c !!!! flnmra = flnmrso !.a extension added by restart_out !!!! flnmrb = flnmrso !.b extension added by restart_out !!!!#endif !!!!c !!!!#if defined (NERSC_VERSION) !!!!c --- This is saved earlier in NERSC_VERSION !!!! continue !!!!#else !!!! call restart_out(nstep,dtime, flnmra,flnmrb, nstep.ge.nstep2) !!!!#endif /* NERSC_VERSION */ !!!!c !!!!c --- set layer thickness (incl.bottom pressure) at u,v points !!!!c --- needed because restart_out may have modified dp !!!!c !!!! call xctilr(dp(1-nbdy,1-nbdy,1,1),1,2*kk, nbdy,nbdy, halo_ps) !!!! call xctilr(dpmixl(1-nbdy,1-nbdy,1),1,2, nbdy,nbdy, halo_ps) !!!!c !!!! margin = nbdy !!!!c !!!! do nm=1,2 !!!!!$OMP PARALLEL DO PRIVATE(j,k,l,i) !!!!!$OMP& SCHEDULE(STATIC,jblk) !!!! do j=1-margin,jj+margin !!!! do l=1,isp(j) !!!! if (nm.eq.mod(nstep+1,2)+1) then !!!! do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) !!!! dpbl(i,j)=dpmixl(i,j,nm) !!!! enddo !!!! endif !!!! do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) !!!! p(i,j,1)=0.0 !!!! do k=1,kk !!!! p(i,j,k+1)=p(i,j,k)+dp(i,j,k,nm) !!!! enddo !!!! enddo !!!! enddo !!!! enddo !!!!!$OMP END PARALLEL DO !!!!c !!!! call dpudpv(dpu(1-nbdy,1-nbdy,1,nm), !!!! & dpv(1-nbdy,1-nbdy,1,nm), !!!! & p,depthu,depthv, max(0,margin-1)) !!!!c !!!! enddo !nm=1,2 !!!!c !!!! call xctmr1(51) !!!! endif ! restrt !!!!c !!!! if (histry .or. hitile .or. hisurf .or. restrt) then !!!! if (mnproc.eq.1) then !!!! write (lp,'(a,i9,a,f9.2,a)') !!!! & ' step',nstep,' day',dtime,' -- archiving completed --' !!!! call flush(lp) !!!! endif !1st tile !!!! endif !!!!c !!!!c --- read next incremental update. !!!!c !!!! if (incflg.ne.0) then !!!! call incupd_rd(dtime) !!!! endif !!!!c !!!! delt1=baclin+baclin !!!!c !!!!#if ! defined(USE_CCSM3) !!!! end_of_run = nstep.ge.nstep2 !!!!#endif !!!!c !!!!c --- at end: output float restart file !!!!c !!!! if (synflt .and. end_of_run) then !!!! call floats_restart !!!! endif !synflt+end_of_run end subroutine HYCOM_Run subroutine HYCOM_Final c c --- end of the run. if (mnproc.eq.1) then print*, 'TW (mod_hycom_wavesonly.F): waves-only run finished' write(nod,'(a)') 'normal stop' call flush(nod) endif !1st tile call xcstop('(normal)') !calls xctmrp stop '(normal)' !won't get here end subroutine HYCOM_Final end module mod_hycom_wavesonly c c c> Revision history: c> c> May 1997 - removed statement "theta(1)=-thbase" after loop 14 c> June 1997 - added loop 60 to fix bug in vertical summation of -diaflx- c> Oct. 1999 - option for krt or kpp mixed layer model - convec and diapfl c> not called for kpp mixing model c> Oct. 1999 - dpbl (boundary layer thickness) is output in addition to c> dpmixl when the kpp mixing model is selected c> May 2000 - conversion to SI units c> Aug. 2000 - added isopycnic (MICOM) vertical coordinate option c> Oct. 2000 - added option for high frequency atmospheric forcing c> Nov. 2000 - archive time stamp is either time step or YYYY_DDD_HH c> Aug. 2002 - added support for multiple tracers c> Nov. 2002 - more basin-wide surface flux statistics c> Dec. 2003 - more basin-wide mean statistics c> Mar. 2005 - more accurate ice statistics c> Jan. 2006 - mod_hycom with HYCOM_Init, HYCOM_Run, HYCOM_Final c> Nov. 2006 - version 2.2 c> Nov. 2006 - added incremental update (for data assimilation) c> Mar. 2007 - added srfhgt