global_cycle  1.12.0
 All Data Structures Files Functions Variables Pages
cycle.f90
Go to the documentation of this file.
1 
4 
101  PROGRAM sfc_drv
102 
103  use mpi
104 
105  IMPLICIT NONE
106 !
107  CHARACTER(LEN=3) :: donst
108  INTEGER :: idim, jdim, lsoil, lugb, iy, im, id, ih, ialb
109  INTEGER :: isot, ivegsrc, lensfc, zsea1_mm, zsea2_mm, ierr
110  INTEGER :: nprocs, myrank, num_threads, num_parthds, max_tasks
111  REAL :: fh, deltsfc, zsea1, zsea2
112  LOGICAL :: use_ufo, do_nsst, do_lndinc, do_sfccycle, frac_grid
113 !
114  namelist/namcyc/ idim,jdim,lsoil,lugb,iy,im,id,ih,fh,&
115  deltsfc,ialb,use_ufo,donst, &
116  do_sfccycle,isot,ivegsrc,zsea1_mm, &
117  zsea2_mm, max_tasks, do_lndinc, frac_grid
118 !
119  DATA idim,jdim,lsoil/96,96,4/
120  DATA iy,im,id,ih,fh/1997,8,2,0,0./
121  DATA lugb/51/, deltsfc/0.0/, ialb/1/, max_tasks/99999/
122  DATA isot/1/, ivegsrc/2/, zsea1_mm/0/, zsea2_mm/0/
123 !
124  CALL mpi_init(ierr)
125  CALL mpi_comm_size(mpi_comm_world, nprocs, ierr)
126  CALL mpi_comm_rank(mpi_comm_world, myrank, ierr)
127 
128  if (myrank==0) call w3tagb('GLOBAL_CYCLE',2018,0179,0055,'NP20')
129 
130  num_threads = num_parthds()
131 
132  print*
133  print*,"STARTING CYCLE PROGRAM ON RANK ", myrank
134  print*,"RUNNING WITH ", nprocs, "TASKS"
135  print*,"AND WITH ", num_threads, " THREADS."
136 
137  use_ufo = .false.
138  donst = "NO"
139  do_lndinc = .false.
140  do_sfccycle = .true.
141  frac_grid = .false.
142 
143  print*
144  print*,"READ NAMCYC NAMELIST."
145 
146  CALL baopenr(36, "fort.36", ierr)
147  READ(36, nml=namcyc)
148 !IF (MYRANK==0) WRITE(6,NAMCYC)
149 
150  IF (max_tasks < 99999 .AND. myrank > (max_tasks - 1)) THEN
151  print*,"USER SPECIFIED MAX NUMBER OF TASKS: ", max_tasks
152  print*,"WILL NOT RUN CYCLE PROGRAM ON RANK: ", myrank
153  goto 333
154  ENDIF
155 
156  lensfc = idim*jdim ! TOTAL NUMBER OF POINTS FOR THE CUBED-SPHERE TILE
157 
158  zsea1 = float(zsea1_mm) / 1000.0 ! CONVERT FROM MM TO METERS
159  zsea2 = float(zsea2_mm) / 1000.0
160 
161  IF (donst == "YES") THEN
162  do_nsst=.true.
163  ELSE
164  do_nsst=.false.
165  ENDIF
166 
167  print*
168  IF (myrank==0) print*,"LUGB,IDIM,JDIM,ISOT,IVEGSRC,LSOIL,DELTSFC,IY,IM,ID,IH,FH: ", &
169  lugb,idim,jdim,isot,ivegsrc,lsoil,deltsfc,iy,im,id,ih,fh
170 
171  CALL sfcdrv(lugb,idim,jdim,lensfc,lsoil,deltsfc, &
172  iy,im,id,ih,fh,ialb, &
173  use_ufo,do_nsst,do_sfccycle,do_lndinc, &
174  frac_grid,zsea1,zsea2,isot,ivegsrc,myrank)
175 
176  print*
177  print*,'CYCLE PROGRAM COMPLETED NORMALLY ON RANK: ', myrank
178 
179  333 CONTINUE
180 
181  CALL mpi_barrier(mpi_comm_world, ierr)
182 
183  if (myrank==0) call w3tage('GLOBAL_CYCLE')
184 
185  CALL mpi_finalize(ierr)
186 
187  stop
188 
189  END PROGRAM sfc_drv
190 
299  SUBROUTINE sfcdrv(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, &
300  iy,im,id,ih,fh,ialb, &
301  use_ufo,do_nsst,do_sfccycle,do_lndinc,&
302  frac_grid,zsea1,zsea2,isot,ivegsrc,myrank)
303 !
304  USE read_write_data
305  use machine
306  USE mpi
307  USE land_increments, ONLY: add_increment_soil, &
308  add_increment_snow, &
309  calculate_landinc_mask, &
310  apply_land_da_adjustments_soil, &
311  apply_land_da_adjustments_snd, &
312  lsm_noah, lsm_noahmp
313 
314  IMPLICIT NONE
315 
316  INTEGER, INTENT(IN) :: idim, jdim, lensfc, lsoil, ialb
317  INTEGER, INTENT(IN) :: lugb, iy, im, id, ih
318  INTEGER, INTENT(IN) :: isot, ivegsrc, myrank
319 
320  LOGICAL, INTENT(IN) :: use_ufo, do_nsst,do_sfccycle
321  LOGICAL, INTENT(IN) :: do_lndinc, frac_grid
322 
323  REAL, INTENT(IN) :: fh, deltsfc, zsea1, zsea2
324 
325  INTEGER, PARAMETER :: nlunit=35
326  INTEGER, PARAMETER :: sz_nml=1
327 
328  CHARACTER(LEN=5) :: tile_num
329  CHARACTER(LEN=500) :: nst_file
330  CHARACTER(LEN=500) :: lnd_soi_file
331  CHARACTER(LEN=4) :: input_nml_file(sz_nml)
332 
333  INTEGER :: i, ierr
334  INTEGER :: i_index(lensfc), j_index(lensfc)
335  INTEGER :: idum(idim,jdim)
336  integer :: num_parthds, num_threads
337 
338  LOGICAL :: is_noahmp
339  INTEGER :: lsm
340 
341  real(kind=kind_io8) :: min_ice(lensfc)
342 
343  REAL :: slmask(lensfc), orog(lensfc)
344  REAL :: sihfcs(lensfc), sicfcs(lensfc)
345  REAL :: sitfcs(lensfc), tsffcs(lensfc)
346  REAL :: swefcs(lensfc), zorfcs(lensfc)
347  REAL :: albfcs(lensfc,4), tg3fcs(lensfc)
348  REAL :: cnpfcs(lensfc), smcfcs(lensfc,lsoil)
349  REAL :: stcfcs(lensfc,lsoil), slifcs(lensfc)
350  REAL :: aisfcs(lensfc), f10m(lensfc)
351  REAL :: vegfcs(lensfc), vetfcs(lensfc)
352  REAL :: sotfcs(lensfc), alffcs(lensfc,2)
353  REAL :: cvfcs(lensfc), cvtfcs(lensfc)
354  REAL :: cvbfcs(lensfc), tprcp(lensfc)
355  REAL :: srflag(lensfc), sndfcs(lensfc)
356  REAL :: slcfcs(lensfc,lsoil), vmxfcs(lensfc)
357  REAL :: vmnfcs(lensfc), t2m(lensfc)
358  REAL :: q2m(lensfc), slpfcs(lensfc)
359  REAL :: absfcs(lensfc), orog_uf(lensfc)
360  REAL :: ustar(lensfc)
361  REAL :: fmm(lensfc), fhh(lensfc)
362  REAL :: rla(lensfc), rlo(lensfc)
363  REAL(KIND=4) :: zsoil(lsoil)
364  REAL :: sig1t(lensfc)
367  REAL, ALLOCATABLE :: stc_bck(:,:), smc_bck(:,:), slc_bck(:,:)
368  REAL, ALLOCATABLE :: slifcs_fg(:), sicfcs_fg(:)
369  INTEGER, ALLOCATABLE :: landinc_mask_fg(:), landinc_mask(:)
370  REAL, ALLOCATABLE :: snd_bck(:), snd_inc(:), swe_bck(:)
371  REAL(KIND=KIND_IO8), ALLOCATABLE :: slmaskl(:), slmaskw(:), landfrac(:)
372 
373  TYPE(nsst_data) :: nsst
374  real, dimension(idim,jdim) :: tf_clm,tf_trd,sal_clm
375  real, dimension(lensfc) :: tf_clm_tile,tf_trd_tile,sal_clm_tile
376  INTEGER :: veg_type_landice
377  INTEGER, DIMENSION(LENSFC) :: stc_updated, slc_updated
378 
379  LOGICAL :: file_exists, do_soi_inc, do_sno_inc
380  CHARACTER(LEN=3) :: rankch
381 
382 !--------------------------------------------------------------------------------
383 ! NST_FILE is the path/name of the gaussian GSI file which contains NSST
384 ! increments.
385 !--------------------------------------------------------------------------------
386 
387  namelist/namsfcd/ nst_file, lnd_soi_file, do_sno_inc
388 
389  DATA nst_file/'NULL'/
390  DATA lnd_soi_file/'NULL'/
391 
392  do_sno_inc = .false.
393  do_soi_inc = .false.
394 
395  sig1t = 0.0 ! Not a dead start!
396 
397  input_nml_file = "NULL"
398 
399  CALL baopenr(37, "fort.37", ierr)
400  READ (37, nml=namsfcd)
401 
402  print*
403  print*,'IN ROUTINE SFCDRV,IDIM=',idim,'JDIM=',jdim,'FH=',fh
404 
405 !--------------------------------------------------------------------------------
406 ! READ THE OROGRAPHY AND GRID POINT LAT/LONS FOR THE CUBED-SPHERE TILE.
407 !--------------------------------------------------------------------------------
408 
409  ALLOCATE(landfrac(lensfc))
410  IF(frac_grid) THEN
411  print*,'- RUNNING WITH FRACTIONAL GRID.'
412  CALL read_lat_lon_orog(rla,rlo,orog,orog_uf,tile_num,idim,jdim,lensfc,landfrac=landfrac)
413  ELSE
414  CALL read_lat_lon_orog(rla,rlo,orog,orog_uf,tile_num,idim,jdim,lensfc)
415  landfrac=-999.9
416  ENDIF
417 
418  DO i = 1, idim
419  idum(i,:) = i
420  ENDDO
421 
422  i_index = reshape(idum, (/lensfc/))
423 
424  DO i = 1, jdim
425  idum(:,i) = i
426  ENDDO
427 
428  j_index = reshape(idum, (/lensfc/) )
429 
430  IF (do_nsst) THEN
431  print*
432  print*,"WILL PROCESS NSST RECORDS."
433  ALLOCATE(nsst%C_0(lensfc))
434  ALLOCATE(nsst%C_D(lensfc))
435  ALLOCATE(nsst%D_CONV(lensfc))
436  ALLOCATE(nsst%DT_COOL(lensfc))
437  ALLOCATE(nsst%IFD(lensfc))
438  ALLOCATE(nsst%QRAIN(lensfc))
439  ALLOCATE(nsst%TREF(lensfc))
440  ALLOCATE(nsst%TFINC(lensfc))
441  ALLOCATE(nsst%W_0(lensfc))
442  ALLOCATE(nsst%W_D(lensfc))
443  ALLOCATE(nsst%XS(lensfc))
444  ALLOCATE(nsst%XT(lensfc))
445  ALLOCATE(nsst%XTTS(lensfc))
446  ALLOCATE(nsst%XU(lensfc))
447  ALLOCATE(nsst%XV(lensfc))
448  ALLOCATE(nsst%XZ(lensfc))
449  ALLOCATE(nsst%XZTS(lensfc))
450  ALLOCATE(nsst%Z_C(lensfc))
451  ALLOCATE(nsst%ZM(lensfc))
452  ALLOCATE(slifcs_fg(lensfc))
453  ALLOCATE(sicfcs_fg(lensfc))
454  ENDIF
455 
456 IF (do_lndinc) THEN
457  IF (frac_grid) THEN
458  print *, 'FATAL ERROR: land increment update does not work with fractional grids.'
459  call mpi_abort(mpi_comm_world, 17, ierr)
460  ENDIF
461  ! identify variables to be updated, and allocate arrays.
462  IF (trim(lnd_soi_file) .NE. "NULL") THEN
463  do_soi_inc = .true.
464  print*
465  print*," APPLYING SOIL INCREMENTS FROM THE GSI"
466  ALLOCATE(stc_bck(lensfc, lsoil), smc_bck(lensfc, lsoil), slc_bck(lensfc,lsoil))
467  ALLOCATE(landinc_mask_fg(lensfc))
468  ENDIF
469  ! FOR NOW, CODE SO CAN DO BOTH, BUT MIGHT NEED TO THINK ABOUT THIS SOME MORE.
470  IF (do_sno_inc) THEN
471  ! ideally, would check here that sfcsub snow DA update is not also requested
472  ! but latter is controlled by fnsol, which is read in within that routine.
473  ! should be done at script level.
474  print*
475  print*," APPLYING SNOW INCREMENTS FROM JEDI"
476  ALLOCATE(snd_bck(lensfc), snd_inc(lensfc), swe_bck(lensfc))
477  ENDIF
478  ! set-up land mask info
479  ALLOCATE(landinc_mask(lensfc))
480  if (ivegsrc == 2) then ! sib
481  veg_type_landice=13
482  else
483  veg_type_landice=15
484  endif
485 ENDIF
486 
487 !--------------------------------------------------------------------------------
488 ! READ THE INPUT SURFACE DATA ON THE CUBED-SPHERE TILE.
489 !--------------------------------------------------------------------------------
490 
491  CALL read_data(lsoil,lensfc,do_nsst,.false.,is_noahmp=is_noahmp, &
492  tsffcs=tsffcs,smcfcs=smcfcs, &
493  swefcs=swefcs,stcfcs=stcfcs,tg3fcs=tg3fcs,zorfcs=zorfcs, &
494  cvfcs=cvfcs, cvbfcs=cvbfcs,cvtfcs=cvtfcs,albfcs=albfcs, &
495  vegfcs=vegfcs,slifcs=slifcs,cnpfcs=cnpfcs,f10m=f10m , &
496  vetfcs=vetfcs,sotfcs=sotfcs,alffcs=alffcs,ustar=ustar , &
497  fmm=fmm ,fhh=fhh ,sihfcs=sihfcs,sicfcs=sicfcs, &
498  sitfcs=sitfcs,tprcp=tprcp ,srflag=srflag,sndfcs=sndfcs, &
499  vmnfcs=vmnfcs,vmxfcs=vmxfcs,slcfcs=slcfcs,slpfcs=slpfcs, &
500  absfcs=absfcs,t2m=t2m ,q2m=q2m ,slmask=slmask, &
501  zsoil=zsoil, nsst=nsst)
502 
503  IF (frac_grid .AND. .NOT. is_noahmp) THEN
504  print *, 'FATAL ERROR: NOAH lsm update does not work with fractional grids.'
505  call mpi_abort(mpi_comm_world, 18, ierr)
506  ENDIF
507 
508  IF (frac_grid .AND. do_sno_inc) THEN
509  print *, 'FATAL ERROR: Snow increment update does not work with fractional grids.'
510  call mpi_abort(mpi_comm_world, 19, ierr)
511  ENDIF
512 
513  IF (is_noahmp .AND. do_sno_inc) THEN
514  print *, 'FATAL ERROR: Snow increment update does not work with NOAH_MP.'
515  call mpi_abort(mpi_comm_world, 29, ierr)
516  ENDIF
517 
518  IF (is_noahmp) THEN
519  lsm=lsm_noahmp
520  ELSE
521  lsm=lsm_noah
522  ENDIF
523 
524  IF (use_ufo) THEN
525  print*
526  print*,'USE UNFILTERED OROGRAPHY.'
527  ELSE
528  orog_uf = 0.0
529  ENDIF
530 
531  DO i=1,lensfc
532  aisfcs(i) = 0.
533  IF(nint(slifcs(i)).EQ.2) aisfcs(i) = 1.
534  ENDDO
535 
536  IF (do_nsst) THEN
537  sicfcs_fg=sicfcs
538  IF (.NOT. do_sfccycle ) THEN
539  print*
540  print*,"FIRST GUESS MASK ADJUSTED BY IFD RECORD"
541  slifcs_fg = slifcs
542  WHERE(nint(nsst%IFD) == 3) slifcs_fg = 2.0
543  ELSE
544  print*
545  print*,"SAVE FIRST GUESS MASK"
546  slifcs_fg = slifcs
547  ENDIF
548  ENDIF
549 
550  ! CALCULATE MASK FOR LAND INCREMENTS
551  IF (do_lndinc) &
552  CALL calculate_landinc_mask(smcfcs(:,1),swefcs, vetfcs, &
553  lensfc,veg_type_landice, landinc_mask)
554 
555 !--------------------------------------------------------------------------------
556 ! UPDATE SURFACE FIELDS.
557 !
558 ! FIRST, SET WATER AND LAND MASKS - SLMASKW/SLMASKL. FOR UNCOUPLED
559 ! (NON-FRACTIONAL) MODE, THESE ARE IDENTICAL TO THE CURRENT
560 ! MASK - '0' WATER; '1' LAND.
561 !--------------------------------------------------------------------------------
562 
563  IF (do_sfccycle) THEN
564  ALLOCATE(slmaskl(lensfc), slmaskw(lensfc))
565 
566  set_mask : IF (frac_grid) THEN
567 
568  DO i=1,lensfc
569  IF(landfrac(i) > 0.0_kind_io8) THEN
570  slmaskl(i) = ceiling(landfrac(i)-1.0e-6_kind_io8)
571  slmaskw(i) = floor(landfrac(i)+1.0e-6_kind_io8)
572  ELSE
573  IF(nint(slmask(i)) == 1) THEN ! If landfrac is zero, this should not happen.
574  ! So, stop processing.
575  print*, 'FATAL ERROR: LAND FRAC AND SLMASK MISMATCH.'
576  CALL mpi_abort(mpi_comm_world, 27, ierr)
577  ELSE
578  slmaskl(i) = 0.0_kind_io8
579  slmaskw(i) = 0.0_kind_io8
580  ENDIF
581  ENDIF
582 
583  ENDDO
584 
585  ELSE
586 
587 ! For running uncoupled (non-fractional grid).
588 
589  DO i=1,lensfc
590  IF(nint(slmask(i)) == 1) THEN
591  slmaskl(i) = 1.0_kind_io8
592  slmaskw(i) = 1.0_kind_io8
593  ELSE
594  slmaskl(i) = 0.0_kind_io8
595  slmaskw(i) = 0.0_kind_io8
596  ENDIF
597  ENDDO
598 
599  ENDIF set_mask
600 
601  DO i=1,lensfc
602  if(nint(slmask(i)) == 0) then
603  min_ice(i) = 0.15_kind_io8
604  else
605  min_ice(i) = 0.0_kind_io8
606  endif
607  ENDDO
608 
609  num_threads = num_parthds()
610  print*
611  print*,"CALL SFCCYCLE TO UPDATE SURFACE FIELDS."
612  CALL sfccycle(lugb,lensfc,lsoil,sig1t,deltsfc, &
613  iy,im,id,ih,fh,rla,rlo, &
614  slmaskl,slmaskw, orog, orog_uf, use_ufo, do_nsst, &
615  sihfcs,sicfcs,sitfcs,sndfcs,slcfcs, &
616  vmnfcs,vmxfcs,slpfcs,absfcs, &
617  tsffcs,swefcs,zorfcs,albfcs,tg3fcs, &
618  cnpfcs,smcfcs,stcfcs,slifcs,aisfcs, &
619  vegfcs,vetfcs,sotfcs,alffcs, &
620  cvfcs,cvbfcs,cvtfcs,myrank,num_threads, nlunit, &
621  sz_nml, input_nml_file, &
622  min_ice, &
623  ialb,isot,ivegsrc,tile_num,i_index,j_index)
624 
625  DEALLOCATE(slmaskl, slmaskw)
626  ENDIF
627 
628 !--------------------------------------------------------------------------------
629 ! IF RUNNING WITH NSST, READ IN GSI FILE WITH THE UPDATED INCREMENTS (ON THE
630 ! GAUSSIAN GRID), INTERPOLATE INCREMENTS TO THE CUBED-SPHERE TILE, AND PERFORM
631 ! REQUIRED ADJUSTMENTS AND QC.
632 !--------------------------------------------------------------------------------
633 
634  IF (do_nsst) THEN
635  IF (nst_file == "NULL") THEN
636  print*
637  print*,"NO GSI FILE. ADJUST IFD FOR FORMER ICE POINTS."
638  DO i = 1, lensfc
639  IF (sicfcs_fg(i) > 0.0 .AND. sicfcs(i) == 0) THEN
640  nsst%IFD(i) = 3.0
641  ENDIF
642  ENDDO
643  nsst%TFINC = 0.0
644  ELSE
645  print*
646  print*,"ADJUST TREF FROM GSI INCREMENT"
647 !
648 ! Get tf climatology at the time
649 !
650  call get_tf_clm(rla,rlo,jdim,idim,iy,im,id,ih,tf_clm,tf_trd)
651  tf_clm_tile(:) = reshape(tf_clm, (/lensfc/) )
652  tf_trd_tile(:) = reshape(tf_trd, (/lensfc/) )
653 !
654 ! Get salinity climatology at the time
655 !
656  call get_sal_clm(rla,rlo,jdim,idim,iy,im,id,ih,sal_clm)
657  sal_clm_tile(:) = reshape(sal_clm, (/lensfc/) )
658 !
659 ! read tf analysis increment generated by GSI
660 !
661  CALL read_gsi_data(nst_file, 'NST')
662 !
663 ! update foundation & surface temperature for NSST
664 !
665  CALL adjust_nsst(rla,rlo,slifcs,slifcs_fg,tsffcs,sitfcs,sicfcs,sicfcs_fg,&
666  stcfcs,nsst,lensfc,lsoil,idim,jdim,zsea1,zsea2, &
667  tf_clm_tile,tf_trd_tile,sal_clm_tile,landfrac,frac_grid)
668  ENDIF
669  ENDIF
670 
671 !--------------------------------------------------------------------------------
672 ! READ IN AND APPLY LAND INCREMENTS FROM THE GSI
673 !--------------------------------------------------------------------------------
674 
675  IF (do_lndinc) THEN
676 
677  ! SNOW INCREMENTS
678  ! do snow first, as temperature updates will use snow analaysis
679  IF (do_sno_inc) THEN
680  ! updates are made to snow depth only over land (and not-land ice).
681  ! SWE is then updated from the snow depth analysis, using the model
682  ! forecast density
683 
684  !--------------------------------------------------------------------------------
685  ! read increments in
686  !--------------------------------------------------------------------------------
687 
688  ! Only coded for DA on native model grid (would always be the case for cycling DA)
689  CALL read_data(lsoil,lensfc,.false.,.true.,sndfcs=snd_inc)
690 
691  !--------------------------------------------------------------------------------
692  ! add increments to state vars
693  !--------------------------------------------------------------------------------
694 
695  ! store background states
696  snd_bck = sndfcs
697  swe_bck = swefcs
698 
699  CALL add_increment_snow(snd_inc,landinc_mask,lensfc,sndfcs)
700 
701  !--------------------------------------------------------------------------------
702  ! make any necessary adjustments to dependent variables
703  !--------------------------------------------------------------------------------
704 
705  CALL apply_land_da_adjustments_snd(lsm, lensfc, landinc_mask, swe_bck, snd_bck, &
706  sndfcs, swefcs)
707 
708  ENDIF
709 
710  ! SOIL INCREMENTS
711  IF (do_soi_inc) THEN
712 
713  !--------------------------------------------------------------------------------
714  ! re-calculate soilsnow mask if snow has been updated.
715  !--------------------------------------------------------------------------------
716 
717  landinc_mask_fg = landinc_mask
718 
719  IF (do_sfccycle .OR. do_sno_inc) THEN
720  CALL calculate_landinc_mask(smcfcs(:,1),swefcs, vetfcs, lensfc, &
721  veg_type_landice, landinc_mask )
722  ENDIF
723 
724  !--------------------------------------------------------------------------------
725  ! read increments in
726  !--------------------------------------------------------------------------------
727 
728  WRITE(rankch, '(I3.3)') (myrank+1)
729 
730  lnd_soi_file = trim(lnd_soi_file) // "." // rankch
731 
732  INQUIRE(file=trim(lnd_soi_file), exist=file_exists)
733  IF (.not. file_exists) then
734  print *, 'FATAL ERROR: land increment update requested, but file does not exist: ', &
735  trim(lnd_soi_file)
736  call mpi_abort(mpi_comm_world, 10, ierr)
737  ENDIF
738 
739  CALL read_gsi_data(lnd_soi_file, 'LND', lsoil=lsoil)
740 
741  !--------------------------------------------------------------------------------
742  ! add increments to state vars
743  !--------------------------------------------------------------------------------
744  ! when applying increments, will often need to adjust other land states in response
745  ! to the changes made. Need to store bacground, apply the increments, then make
746  ! secondart adjustments. When updating more than one state, be careful about the
747  ! order if increments and secondary adjustments.
748 
749  ! store background states
750  stc_bck = stcfcs
751  smc_bck = smcfcs
752  slc_bck = slcfcs
753 
754  ! below updates [STC/SMC/STC]FCS to hold the analysis
755  CALL add_increment_soil(rla,rlo,stcfcs,smcfcs,slcfcs,stc_updated, slc_updated, &
756  landinc_mask,landinc_mask_fg,lensfc,lsoil,idim,jdim,lsm,myrank)
757 
758  !--------------------------------------------------------------------------------
759  ! make any necessary adjustments to dependent variables
760  !--------------------------------------------------------------------------------
761 
762 
763  CALL apply_land_da_adjustments_soil(lsm, isot, ivegsrc,lensfc, lsoil, &
764  sotfcs, landinc_mask_fg, stc_bck, stcfcs, smcfcs, slcfcs, stc_updated, &
765  slc_updated,zsoil)
766 
767  ENDIF ! soil increments
768 
769 !--------------------------------------------------------------------------------
770 ! clean up
771 !--------------------------------------------------------------------------------
772 
773  ! to do - save and write out STC_INC? (soil temperature increments)
774  IF(ALLOCATED(landinc_mask_fg)) DEALLOCATE(landinc_mask_fg)
775  IF(ALLOCATED(landinc_mask)) DEALLOCATE(landinc_mask)
776  IF(ALLOCATED(stc_bck)) DEALLOCATE(stc_bck)
777  IF(ALLOCATED(smc_bck)) DEALLOCATE(smc_bck)
778  IF(ALLOCATED(slc_bck)) DEALLOCATE(slc_bck)
779  IF(ALLOCATED(snd_bck)) DEALLOCATE(snd_bck)
780  IF(ALLOCATED(swe_bck)) DEALLOCATE(swe_bck)
781  IF(ALLOCATED(snd_inc)) DEALLOCATE(snd_inc)
782 
783  ENDIF
784 !--------------------------------------------------------------------------------
785 ! WRITE OUT UPDATED SURFACE DATA ON THE CUBED-SPHERE TILE.
786 !--------------------------------------------------------------------------------
787 
788  IF (lsm==lsm_noahmp) THEN
789 
790  CALL write_data(lensfc,idim,jdim,lsoil,do_nsst,nsst,vegfcs=vegfcs, &
791  slcfcs=slcfcs,smcfcs=smcfcs,stcfcs=stcfcs,&
792  sicfcs=sicfcs,sihfcs=sihfcs)
793 
794  ELSEIF (lsm==lsm_noah) THEN
795 
796  CALL write_data(lensfc,idim,jdim,lsoil, &
797  do_nsst,nsst,slifcs=slifcs,tsffcs=tsffcs,vegfcs=vegfcs, &
798  swefcs=swefcs,tg3fcs=tg3fcs,zorfcs=zorfcs, &
799  albfcs=albfcs,alffcs=alffcs,cnpfcs=cnpfcs, &
800  f10m=f10m,t2m=t2m,q2m=q2m,vetfcs=vetfcs, &
801  sotfcs=sotfcs,ustar=ustar,fmm=fmm,fhh=fhh, &
802  sicfcs=sicfcs,sihfcs=sihfcs,sitfcs=sitfcs,tprcp=tprcp, &
803  srflag=srflag,swdfcs=sndfcs,vmnfcs=vmnfcs, &
804  vmxfcs=vmxfcs,slpfcs=slpfcs,absfcs=absfcs, &
805  slcfcs=slcfcs,smcfcs=smcfcs,stcfcs=stcfcs)
806 
807  ENDIF
808 
809  IF (do_nsst) THEN
810  DEALLOCATE(nsst%C_0)
811  DEALLOCATE(nsst%C_D)
812  DEALLOCATE(nsst%D_CONV)
813  DEALLOCATE(nsst%DT_COOL)
814  DEALLOCATE(nsst%IFD)
815  DEALLOCATE(nsst%QRAIN)
816  DEALLOCATE(nsst%TREF)
817  DEALLOCATE(nsst%TFINC)
818  DEALLOCATE(nsst%W_0)
819  DEALLOCATE(nsst%W_D)
820  DEALLOCATE(nsst%XS)
821  DEALLOCATE(nsst%XT)
822  DEALLOCATE(nsst%XTTS)
823  DEALLOCATE(nsst%XU)
824  DEALLOCATE(nsst%XV)
825  DEALLOCATE(nsst%XZ)
826  DEALLOCATE(nsst%XZTS)
827  DEALLOCATE(nsst%Z_C)
828  DEALLOCATE(nsst%ZM)
829  DEALLOCATE(slifcs_fg)
830  DEALLOCATE(sicfcs_fg)
831  ENDIF
832 
833  RETURN
834 
835  END SUBROUTINE sfcdrv
836 
868  SUBROUTINE adjust_nsst(RLA,RLO,SLMSK_TILE,SLMSK_FG_TILE,SKINT_TILE,&
869  sicet_tile,sice_tile,sice_fg_tile,soilt_tile,nsst, &
870  lensfc,lsoil,idim,jdim,zsea1,zsea2, &
871  tf_clm_tile,tf_trd_tile,sal_clm_tile,landfrac, &
872  frac_grid)
873 
874  USE utils
875  USE gdswzd_mod
876  USE read_write_data, ONLY : idim_gaus, jdim_gaus, &
877  slmsk_gaus, dtref_gaus, &
878  nsst_data
879 
880  USE mpi
881 
882  IMPLICIT NONE
883 
884  INTEGER, INTENT(IN) :: lensfc, lsoil, idim, jdim
885 
886  LOGICAL, INTENT(IN) :: frac_grid
887 
888  REAL, INTENT(IN) :: slmsk_tile(lensfc), slmsk_fg_tile(lensfc), landfrac(lensfc)
889  real, intent(in) :: tf_clm_tile(lensfc),tf_trd_tile(lensfc),sal_clm_tile(lensfc)
890  REAL, INTENT(IN) :: zsea1, zsea2,sice_tile(lensfc),sice_fg_tile(lensfc)
891  REAL, INTENT(IN) :: rla(lensfc), rlo(lensfc)
892  REAL, INTENT(INOUT) :: skint_tile(lensfc)
893  REAL, INTENT(INOUT) :: sicet_tile(lensfc),soilt_tile(lensfc,lsoil)
894 
895  TYPE(nsst_data) :: nsst
896 
897  REAL, PARAMETER :: tmax=313.0,tzero=273.16
898 
899  INTEGER :: iopt, nret, kgds_gaus(200)
900  INTEGER :: igaus, jgaus, ij, ii, jj, iii, jjj, krad
901  INTEGER :: istart, iend, jstart, jend
902 !INTEGER :: MASK_TILE, MASK_FG_TILE
903  INTEGER,allocatable :: mask_tile(:),mask_fg_tile(:)
904  INTEGER :: itile, jtile
905  INTEGER :: max_search, j, ierr
906  INTEGER :: igausp1, jgausp1
907  integer :: nintp,nsearched,nice,nland
908  integer :: nfill,nfill_tice,nfill_clm
909  integer :: nset_thaw,nset_thaw_s,nset_thaw_i,nset_thaw_c
910 
911  INTEGER, ALLOCATABLE :: id1(:,:), id2(:,:), jdc(:,:)
912 
913  LOGICAL :: is_ice
914 
915  real :: tfreez
916  REAL :: tref_save,wsum,tf_ice,tf_thaw
917  REAL :: fill, dtzm, gaus_res_km, dtref
918  REAL, ALLOCATABLE :: xpts(:), ypts(:), lats(:), lons(:)
919  REAL, ALLOCATABLE :: dum2d(:,:), lats_rad(:), lons_rad(:)
920  REAL, ALLOCATABLE :: agrid(:,:,:), s2c(:,:,:)
921 
922  kgds_gaus = 0
923  kgds_gaus(1) = 4 ! OCT 6 - TYPE OF GRID (GAUSSIAN)
924  kgds_gaus(2) = idim_gaus ! OCT 7-8 - # PTS ON LATITUDE CIRCLE
925  kgds_gaus(3) = jdim_gaus
926  kgds_gaus(4) = 90000 ! OCT 11-13 - LAT OF ORIGIN
927  kgds_gaus(5) = 0 ! OCT 14-16 - LON OF ORIGIN
928  kgds_gaus(6) = 128 ! OCT 17 - RESOLUTION FLAG
929  kgds_gaus(7) = -90000 ! OCT 18-20 - LAT OF EXTREME POINT
930  kgds_gaus(8) = nint(-360000./float(idim_gaus)) ! OCT 21-23 - LON OF EXTREME POINT
931  kgds_gaus(9) = nint((360.0 / float(idim_gaus))*1000.0)
932  ! OCT 24-25 - LONGITUDE DIRECTION INCR.
933  kgds_gaus(10) = jdim_gaus/2 ! OCT 26-27 - NUMBER OF CIRCLES POLE TO EQUATOR
934  kgds_gaus(12) = 255 ! OCT 29 - RESERVED
935  kgds_gaus(20) = 255 ! OCT 5 - NOT USED, SET TO 255
936 
937  print*
938  print*,'ADJUST NSST USING GSI INCREMENTS ON GAUSSIAN GRID'
939 
940 !----------------------------------------------------------------------
941 ! CALL GDSWZD TO COMPUTE THE LAT/LON OF EACH GSI GAUSSIAN GRID POINT.
942 !----------------------------------------------------------------------
943 
944  iopt = 0
945  fill = -9999.
946  ALLOCATE(xpts(idim_gaus*jdim_gaus))
947  ALLOCATE(ypts(idim_gaus*jdim_gaus))
948  ALLOCATE(lats(idim_gaus*jdim_gaus))
949  ALLOCATE(lons(idim_gaus*jdim_gaus))
950  xpts = fill
951  ypts = fill
952  lats = fill
953  lons = fill
954 
955  CALL gdswzd(kgds_gaus,iopt,(idim_gaus*jdim_gaus),fill,xpts,ypts,lons,lats,nret)
956 
957  IF (nret /= (idim_gaus*jdim_gaus)) THEN
958  print*,'FATAL ERROR: PROBLEM IN GDSWZD. STOP.'
959  CALL mpi_abort(mpi_comm_world, 12, ierr)
960  ENDIF
961 
962  DEALLOCATE (xpts, ypts)
963 
964  ALLOCATE(dum2d(idim_gaus,jdim_gaus))
965  dum2d = reshape(lats, (/idim_gaus,jdim_gaus/) )
966  DEALLOCATE(lats)
967 
968  ALLOCATE(lats_rad(jdim_gaus))
969  DO j = 1, jdim_gaus
970  lats_rad(j) = dum2d(1,jdim_gaus-j+1) * 3.1415926 / 180.0
971  ENDDO
972 
973  dum2d = reshape(lons, (/idim_gaus,jdim_gaus/) )
974  DEALLOCATE(lons)
975  ALLOCATE(lons_rad(idim_gaus))
976  lons_rad = dum2d(:,1) * 3.1415926 / 180.0
977 
978  DEALLOCATE(dum2d)
979 
980  ALLOCATE(agrid(idim,jdim,2))
981  agrid(:,:,1) = reshape(rlo, (/idim,jdim/) )
982  agrid(:,:,2) = reshape(rla, (/idim,jdim/) )
983  agrid = agrid * 3.1415926 / 180.0
984 
985  ALLOCATE(id1(idim,jdim))
986  ALLOCATE(id2(idim,jdim))
987  ALLOCATE(jdc(idim,jdim))
988  ALLOCATE(s2c(idim,jdim,4))
989 
990 !----------------------------------------------------------------------
991 ! COMPUTE BILINEAR WEIGHTS FOR EACH MODEL POINT FROM THE NEAREST
992 ! FOUR GSI/GAUSSIAN POINTS. DOES NOT ACCOUNT FOR MASK. THAT
993 ! HAPPENS LATER.
994 !----------------------------------------------------------------------
995 
996  CALL remap_coef( 1, idim, 1, jdim, idim_gaus, jdim_gaus, &
997  lons_rad, lats_rad, id1, id2, jdc, s2c, agrid )
998 
999  DEALLOCATE(lons_rad, lats_rad, agrid)
1000 
1001 !----------------------------------------------------------------------
1002 ! THE MAXIMUM DISTANCE TO SEARCH IS 500 KM. HOW MANY GAUSSIAN
1003 ! GRID LENGTHS IS THAT?
1004 !----------------------------------------------------------------------
1005 
1006  gaus_res_km = 360.0 / idim_gaus * 111.0
1007  max_search = ceiling(500.0/gaus_res_km)
1008 
1009  print*
1010  print*,'MAXIMUM SEARCH IS ',max_search, ' GAUSSIAN POINTS.'
1011  print*
1012 
1013 !
1014 ! Initialize variables for counts statitics to be zeros
1015 !
1016  nintp = 0
1017  nset_thaw = 0
1018  nset_thaw_s = 0
1019  nset_thaw_i = 0
1020  nset_thaw_c = 0
1021  nsearched = 0
1022  nfill = 0
1023  nfill_tice = 0
1024  nfill_clm = 0
1025  nice = 0
1026  nland = 0
1027 !----------------------------------------------------------------------
1028 ! TREF INCREMENT WILL BE OUTPUT. INITIALIZE TO ZERO.
1029 !----------------------------------------------------------------------
1030 
1031  nsst%TFINC = 0.0
1032 
1033  allocate(mask_tile(lensfc))
1034  allocate(mask_fg_tile(lensfc))
1035 
1036  IF(.NOT. frac_grid) THEN
1037  mask_tile = nint(slmsk_tile)
1038  mask_fg_tile = nint(slmsk_fg_tile)
1039  ELSE
1040  mask_tile=0
1041  WHERE(sice_tile > 0.0) mask_tile=2
1042  WHERE(landfrac == 1.0) mask_tile=1
1043  mask_fg_tile=0
1044  WHERE(sice_fg_tile > 0.0) mask_fg_tile=2
1045  WHERE(landfrac == 1.0) mask_fg_tile=1
1046  ENDIF
1047 
1048  ij_loop : DO ij = 1, lensfc
1049 
1050 !
1051 ! when sea ice exists, get salinity dependent water temperature
1052 !
1053  tf_ice = tfreez(sal_clm_tile(ij)) + tzero
1054 !----------------------------------------------------------------------
1055 ! SKIP LAND POINTS. NSST NOT APPLIED AT LAND.
1056 !----------------------------------------------------------------------
1057 
1058  IF (mask_tile(ij) == 1) THEN
1059  nland = nland + 1
1060  cycle ij_loop
1061  ENDIF
1062 
1063 !
1064 ! these are ice points. set tref to tf_ice and update tmpsfc.
1065 !
1066  if (mask_tile(ij) == 2) then
1067  nsst%tref(ij)=tf_ice ! water part tmp set
1068  skint_tile(ij)=(1.0-sice_tile(ij))*nsst%tref(ij)+sice_tile(ij)*sicet_tile(ij)
1069  nice = nice + 1
1070  cycle ij_loop
1071  endif
1072 
1073 !
1074 ! Get i,j index on array of (idim,jdim) from known ij
1075 !
1076  jtile = (ij-1) / idim + 1
1077  itile = mod(ij,idim)
1078  IF (itile==0) itile = idim
1079 
1080 !----------------------------------------------------------------------
1081 ! IF THE MODEL POINT WAS ICE COVERED, BUT IS NOW OPEN WATER, SET
1082 ! TREF TO searched adjascent open water onea, if failed the search, set to
1083 ! weighted average of tf_ice and tf_clm. For NSST vars, set xz TO '30' AND ALL OTHER FIELDS TO ZERO.
1084 !----------------------------------------------------------------------
1085 
1086  IF (mask_fg_tile(ij) == 2 .AND. mask_tile(ij) == 0) THEN
1087 !
1088 ! set background for the thaw (just melted water) situation
1089 !
1090  call tf_thaw_set(nsst%tref,mask_fg_tile,itile,jtile,tf_ice,tf_clm_tile(ij),tf_thaw,idim,jdim, &
1091  nset_thaw_s,nset_thaw_i,nset_thaw_c)
1092  call nsst_water_reset(nsst,ij,tf_thaw)
1093  nset_thaw = nset_thaw + 1
1094  ENDIF
1095 
1096 !----------------------------------------------------------------------
1097 ! THESE ARE POINTS THAT ARE OPEN WATER AND WERE OPEN WATER PRIOR
1098 ! TO ANY ICE UPDATE BY SFCCYCLE. UPDATE TREF AND SKIN TEMP.
1099 ! AT OPEN WATER POINTS, THE SEA ICE TEMPERATURE (SICET_TILE) AND
1100 ! SOIL COLUMN TEMPERATURE (SOILT_TILE) ARE SET TO THE SKIN TEMP.
1101 ! IT IS SIMPLY A FILLER VALUE. THESE FIELDS ARE NOT USED AT
1102 ! OPEN WATER POINTS.
1103 !----------------------------------------------------------------------
1104 !----------------------------------------------------------------------
1105 ! SEE IF ANY OF THE NEAREST GSI POINTS MASK AREA OPEN WATER.
1106 ! IF SO, APPLY NSST INCREMENT USING BILINEAR INTERPOLATION.
1107 !----------------------------------------------------------------------
1108 
1109  igaus = id1(itile,jtile)
1110  jgaus = jdc(itile,jtile)
1111  igausp1 = id2(itile,jtile)
1112  jgausp1 = jdc(itile,jtile)+1
1113 
1114  IF (slmsk_gaus(igaus,jgaus) == 0 .OR. &
1115  slmsk_gaus(igausp1,jgaus) == 0 .OR. &
1116  slmsk_gaus(igausp1,jgausp1) == 0 .OR. &
1117  slmsk_gaus(igaus,jgausp1) == 0) THEN
1118 
1119  dtref = 0.0
1120  wsum = 0.0
1121 
1122  IF (slmsk_gaus(igaus,jgaus) == 0) THEN
1123  dtref = dtref + (s2c(itile,jtile,1) * dtref_gaus(igaus,jgaus))
1124  wsum = wsum + s2c(itile,jtile,1)
1125  ENDIF
1126 
1127  IF (slmsk_gaus(igausp1,jgaus) == 0) THEN
1128  dtref = dtref + (s2c(itile,jtile,2) * dtref_gaus(igausp1,jgaus))
1129  wsum = wsum + s2c(itile,jtile,2)
1130  ENDIF
1131 
1132  IF (slmsk_gaus(igausp1,jgausp1) == 0) THEN
1133  dtref = dtref + (s2c(itile,jtile,3) * dtref_gaus(igausp1,jgausp1))
1134  wsum = wsum + s2c(itile,jtile,3)
1135  ENDIF
1136 
1137  IF (slmsk_gaus(igaus,jgausp1) == 0) THEN
1138  dtref = dtref + (s2c(itile,jtile,4) * dtref_gaus(igaus,jgausp1))
1139  wsum = wsum + s2c(itile,jtile,4)
1140  ENDIF
1141 
1142  nintp = nintp + 1
1143  dtref = dtref / wsum
1144 
1145  tref_save = nsst%TREF(ij)
1146  nsst%TREF(ij) = nsst%TREF(ij) + dtref
1147  nsst%TREF(ij) = max(nsst%TREF(ij), tf_ice)
1148  nsst%TREF(ij) = min(nsst%TREF(ij), tmax)
1149  nsst%TFINC(ij) = nsst%TREF(ij) - tref_save
1150 
1151  CALL dtzm_point(nsst%XT(ij),nsst%XZ(ij),nsst%DT_COOL(ij), &
1152  nsst%Z_C(ij),zsea1,zsea2,dtzm)
1153 
1154  skint_tile(ij) = nsst%TREF(ij) + dtzm
1155  skint_tile(ij) = max(skint_tile(ij), tf_ice)
1156  skint_tile(ij) = min(skint_tile(ij), tmax)
1157 
1158  sicet_tile(ij) = skint_tile(ij)
1159 ! Under fractional grids, soilt is used at points with at
1160 ! least some land.
1161  IF(.NOT. frac_grid) soilt_tile(ij,:) = skint_tile(ij)
1162 
1163 !----------------------------------------------------------------------
1164 ! NO NEARBY GSI/GAUSSIAN OPEN WATER POINTS. PERFORM A SPIRAL SEARCH TO
1165 ! FIND NEAREST NON-LAND POINT ON GSI/GAUSSIAN GRID.
1166 !----------------------------------------------------------------------
1167 
1168  ELSE
1169 
1170  is_ice = .false.
1171 
1172  DO krad = 1, max_search
1173 
1174  istart = igaus - krad
1175  iend = igaus + krad
1176  jstart = jgaus - krad
1177  jend = jgaus + krad
1178 
1179  DO jj = jstart, jend
1180  DO ii = istart, iend
1181 
1182  IF((jj == jstart) .OR. (jj == jend) .OR. &
1183  (ii == istart) .OR. (ii == iend)) THEN
1184 
1185  IF ((jj >= 1) .AND. (jj <= jdim_gaus)) THEN
1186 
1187  jjj = jj
1188  IF (ii <= 0) THEN
1189  iii = idim_gaus + ii
1190  ELSE IF (ii >= (idim_gaus+1)) THEN
1191  iii = ii - idim_gaus
1192  ELSE
1193  iii = ii
1194  END IF
1195 
1196 !----------------------------------------------------------------------
1197 ! SEE IF NEARBY POINTS ARE SEA ICE. IF THEY ARE, AND THE SEARCH FOR
1198 ! A GAUSSIAN GRID OPEN WATER POINT FAILS, THEN TREF WILL BE SET TO
1199 ! FREEZING BELOW.
1200 !----------------------------------------------------------------------
1201 
1202  IF (krad <= 2 .AND. slmsk_gaus(iii,jjj) == 2) is_ice = .true.
1203 
1204  IF (slmsk_gaus(iii,jjj) == 0) THEN
1205 
1206 ! PRINT*,'MISMATCH AT TILE POINT ',ITILE,JTILE
1207 ! PRINT*,'UPDATE TREF USING GSI INCREMENT AT ',III,JJJ,DTREF_GAUS(III,JJJ)
1208  nsearched = nsearched + 1
1209 
1210  tref_save = nsst%TREF(ij)
1211  nsst%TREF(ij ) = nsst%TREF(ij) + dtref_gaus(iii,jjj)
1212  nsst%TREF(ij) = max(nsst%TREF(ij), tf_ice)
1213  nsst%TREF(ij) = min(nsst%TREF(ij), tmax)
1214  nsst%TFINC(ij) = nsst%TREF(ij) - tref_save
1215 
1216  CALL dtzm_point(nsst%XT(ij),nsst%XZ(ij),nsst%DT_COOL(ij), &
1217  nsst%Z_C(ij),zsea1,zsea2,dtzm)
1218 
1219  skint_tile(ij) = nsst%TREF(ij) + dtzm
1220  skint_tile(ij) = max(skint_tile(ij), tf_ice)
1221  skint_tile(ij) = min(skint_tile(ij), tmax)
1222 
1223  sicet_tile(ij) = skint_tile(ij)
1224  IF(.NOT. frac_grid) soilt_tile(ij,:) = skint_tile(ij)
1225  cycle ij_loop
1226 
1227  ENDIF ! GSI/Gaussian mask is open water
1228 
1229  ENDIF
1230 
1231  ENDIF
1232 
1233  ENDDO
1234  ENDDO
1235 
1236  ENDDO ! KRAD LOOP
1237 
1238 !----------------------------------------------------------------------
1239 ! THE SEARCH FAILED. IF THERE IS NEARBY ICE, SET TREF TO FREEZING.
1240 ! ELSE UPDATE TREF BASED ON THE ANNUAL SST CYCLE.
1241 !----------------------------------------------------------------------
1242 
1243 ! PRINT*,'WARNING !!!!!! SEARCH FAILED AT TILE POINT ',ITILE,JTILE
1244 
1245  nfill = nfill + 1
1246  IF (is_ice) THEN
1247  nsst%TREF(ij) = tf_ice
1248 ! PRINT*,"NEARBY ICE. SET TREF TO FREEZING"
1249  nfill_tice = nfill_tice + 1
1250  ELSE
1251  tref_save = nsst%TREF(ij)
1252  nsst%TREF(ij) = nsst%TREF(ij) + tf_trd_tile(ij)
1253  nsst%TREF(ij) = max(nsst%TREF(ij), tf_ice)
1254  nsst%TREF(ij) = min(nsst%TREF(ij), tmax)
1255  nsst%TFINC(ij) = nsst%TREF(ij) - tref_save
1256 ! PRINT*,'UPDATE TREF FROM SST CLIMO ',DTREF
1257  nfill_clm = nfill_clm + 1
1258  ENDIF
1259 
1260  CALL dtzm_point(nsst%XT(ij),nsst%XZ(ij),nsst%DT_COOL(ij), &
1261  nsst%Z_C(ij),zsea1,zsea2,dtzm)
1262 
1263  skint_tile(ij) = nsst%TREF(ij) + dtzm
1264  skint_tile(ij) = max(skint_tile(ij), tf_ice)
1265  skint_tile(ij) = min(skint_tile(ij), tmax)
1266 
1267  sicet_tile(ij) = skint_tile(ij)
1268  IF (.NOT. frac_grid) soilt_tile(ij,:) = skint_tile(ij)
1269 
1270  ENDIF ! NEARBY GAUSSIAN POINTS ARE OPEN WATER?
1271 
1272  ENDDO ij_loop
1273 
1274  write(*,'(a)') 'statistics of grids number processed for tile : '
1275  write(*,'(a,I8)') ' nintp = ',nintp
1276  write(*,'(a,4I8)') 'nset_thaw,nset_thaw_s,nset_thaw_i,nset_thaw_c =',nset_thaw,nset_thaw_s,nset_thaw_i,nset_thaw_c
1277  write(*,'(a,I8)') ' nsearched = ',nsearched
1278  write(*,'(a,3I6)') ' nfill,nfill_tice,nfill_clm = ',nfill,nfill_tice,nfill_clm
1279  write(*,'(a,I8)') ' nice = ',nice
1280  write(*,'(a,I8)') ' nland = ',nland
1281 
1282  DEALLOCATE(id1, id2, jdc, s2c, mask_tile, mask_fg_tile)
1283 
1284  END SUBROUTINE adjust_nsst
1285 
1296  SUBROUTINE climo_trend(LATITUDE, MON, DAY, DELTSFC, DTREF)
1297  IMPLICIT NONE
1298 
1299  INTEGER, INTENT(IN) :: mon, day
1300 
1301  REAL, INTENT(IN) :: latitude, deltsfc
1302  REAL, INTENT(OUT) :: dtref
1303 
1304  INTEGER :: num_days(12), mon2, mon1
1305 
1306  REAL, TARGET :: sst_80_90(12)
1307  REAL, TARGET :: sst_70_80(12)
1308  REAL, TARGET :: sst_60_70(12)
1309  REAL, TARGET :: sst_50_60(12)
1310  REAL, TARGET :: sst_40_50(12)
1311  REAL, TARGET :: sst_30_40(12)
1312  REAL, TARGET :: sst_20_30(12)
1313  REAL, TARGET :: sst_10_20(12)
1314  REAL, TARGET :: sst_00_10(12)
1315  REAL, TARGET :: sst_m10_00(12)
1316  REAL, TARGET :: sst_m20_m10(12)
1317  REAL, TARGET :: sst_m30_m20(12)
1318  REAL, TARGET :: sst_m40_m30(12)
1319  REAL, TARGET :: sst_m50_m40(12)
1320  REAL, TARGET :: sst_m60_m50(12)
1321  REAL, TARGET :: sst_m70_m60(12)
1322  REAL, TARGET :: sst_m80_m70(12)
1323  REAL, TARGET :: sst_m90_m80(12)
1324 
1325  REAL, POINTER :: sst(:)
1326 
1327  DATA num_days /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
1328 
1329  DATA sst_80_90 /271.466, 271.458, 271.448, 271.445, 271.519, 271.636, &
1330  272.023, 272.066, 272.001, 271.698, 271.510, 271.472/
1331 
1332  DATA sst_70_80 /272.149, 272.103, 272.095, 272.126, 272.360, 272.988, &
1333  274.061, 274.868, 274.415, 273.201, 272.468, 272.268/
1334 
1335  DATA sst_60_70 /274.240, 274.019, 273.988, 274.185, 275.104, 276.875, &
1336  279.005, 280.172, 279.396, 277.586, 275.818, 274.803/
1337 
1338  DATA sst_50_60 /277.277, 276.935, 277.021, 277.531, 279.100, 281.357, &
1339  283.735, 285.171, 284.399, 282.328, 279.918, 278.199/
1340 
1341  DATA sst_40_50 /281.321, 280.721, 280.850, 281.820, 283.958, 286.588, &
1342  289.195, 290.873, 290.014, 287.652, 284.898, 282.735/
1343 
1344  DATA sst_30_40 /289.189, 288.519, 288.687, 289.648, 291.547, 293.904, &
1345  296.110, 297.319, 296.816, 295.225, 292.908, 290.743/
1346 
1347  DATA sst_20_30 /294.807, 294.348, 294.710, 295.714, 297.224, 298.703, &
1348  299.682, 300.127, 300.099, 299.455, 297.953, 296.177/
1349 
1350  DATA sst_10_20 /298.878, 298.720, 299.033, 299.707, 300.431, 300.709, &
1351  300.814, 300.976, 301.174, 301.145, 300.587, 299.694/
1352 
1353  DATA sst_00_10 /300.415, 300.548, 300.939, 301.365, 301.505, 301.141, &
1354  300.779, 300.660, 300.818, 300.994, 300.941, 300.675/
1355 
1356  DATA sst_m10_00 /300.226, 300.558, 300.914, 301.047, 300.645, 299.870, &
1357  299.114, 298.751, 298.875, 299.294, 299.721, 299.989/
1358 
1359  DATA sst_m20_m10 /299.547, 299.985, 300.056, 299.676, 298.841, 297.788, &
1360  296.893, 296.491, 296.687, 297.355, 298.220, 298.964/
1361 
1362  DATA sst_m30_m20 /297.524, 298.073, 297.897, 297.088, 295.846, 294.520, &
1363  293.525, 293.087, 293.217, 293.951, 295.047, 296.363/
1364 
1365  DATA sst_m40_m30 /293.054, 293.765, 293.468, 292.447, 291.128, 289.781, &
1366  288.773, 288.239, 288.203, 288.794, 289.947, 291.553/
1367 
1368  DATA sst_m50_m40 /285.052, 285.599, 285.426, 284.681, 283.761, 282.826, &
1369  282.138, 281.730, 281.659, 281.965, 282.768, 283.961/
1370 
1371  DATA sst_m60_m50 /277.818, 278.174, 277.991, 277.455, 276.824, 276.229, &
1372  275.817, 275.585, 275.560, 275.687, 276.142, 276.968/
1373 
1374  DATA sst_m70_m60 /273.436, 273.793, 273.451, 272.813, 272.349, 272.048, &
1375  271.901, 271.838, 271.845, 271.889, 272.080, 272.607/
1376 
1377  DATA sst_m80_m70 /271.579, 271.578, 271.471, 271.407, 271.392, 271.391, &
1378  271.390, 271.391, 271.394, 271.401, 271.422, 271.486/
1379 
1380  DATA sst_m90_m80 /271.350, 271.350, 271.350, 271.350, 271.350, 271.350, &
1381  271.350, 271.350, 271.350, 271.350, 271.350, 271.350/
1382 
1383  nullify(sst)
1384  IF (latitude > 80.0) THEN
1385  sst => sst_80_90
1386  ELSEIF (latitude > 70.0) THEN
1387  sst => sst_70_80
1388  ELSEIF (latitude > 60.0) THEN
1389  sst => sst_60_70
1390  ELSEIF (latitude > 50.0) THEN
1391  sst => sst_50_60
1392  ELSEIF (latitude > 40.0) THEN
1393  sst => sst_40_50
1394  ELSEIF (latitude > 30.0) THEN
1395  sst => sst_30_40
1396  ELSEIF (latitude > 20.0) THEN
1397  sst => sst_20_30
1398  ELSEIF (latitude > 10.0) THEN
1399  sst => sst_10_20
1400  ELSEIF (latitude > 0.0) THEN
1401  sst => sst_00_10
1402  ELSEIF (latitude > -10.0) THEN
1403  sst => sst_m10_00
1404  ELSEIF (latitude > -20.0) THEN
1405  sst => sst_m20_m10
1406  ELSEIF (latitude > -30.0) THEN
1407  sst => sst_m30_m20
1408  ELSEIF (latitude > -40.0) THEN
1409  sst => sst_m40_m30
1410  ELSEIF (latitude > -50.0) THEN
1411  sst => sst_m50_m40
1412  ELSEIF (latitude > -60.0) THEN
1413  sst => sst_m60_m50
1414  ELSEIF (latitude > -70.0) THEN
1415  sst => sst_m70_m60
1416  ELSEIF (latitude > -80.0) THEN
1417  sst => sst_m80_m70
1418  ELSE
1419  sst => sst_m90_m80
1420  END IF
1421 
1422  IF (day >= 15) THEN
1423  mon2 = mon+1
1424  IF(mon2 == 13) mon2 = 1
1425  mon1 = mon
1426  dtref = (sst(mon2) - sst(mon1)) / num_days(mon1)
1427  ELSE
1428  mon1 = mon - 1
1429  IF (mon1 == 0) mon1=12
1430  mon2 = mon
1431  dtref = (sst(mon2) - sst(mon1)) / num_days(mon1)
1432  ENDIF
1433 
1434  dtref = dtref * (deltsfc / 24.0)
1435 
1436  END SUBROUTINE climo_trend
1437 
1449  SUBROUTINE dtzm_point(XT,XZ,DT_COOL,ZC,Z1,Z2,DTZM)
1450  implicit none
1451 
1452  real, intent(in) :: xt,xz,dt_cool,zc,z1,z2
1453  real, intent(out) :: dtzm
1454 
1455  real, parameter :: zero = 0.0
1456  real, parameter :: one = 1.0
1457  real, parameter :: half = 0.5
1458  real :: dt_warm,dtw,dtc
1459 !
1460 ! get the mean warming in the range of z=z1 to z=z2
1461 !
1462  dtw = zero
1463  if ( xt > zero ) then
1464  dt_warm = (xt+xt)/xz ! Tw(0)
1465  if ( z1 < z2) then
1466  if ( z2 < xz ) then
1467  dtw = dt_warm*(one-(z1+z2)/(xz+xz))
1468  elseif ( z1 < xz .and. z2 >= xz ) then
1469  dtw = half*(one-z1/xz)*dt_warm*(xz-z1)/(z2-z1)
1470  endif
1471  elseif ( z1 == z2 ) then
1472  if ( z1 < xz ) then
1473  dtw = dt_warm*(one-z1/xz)
1474  endif
1475  endif
1476  endif
1477 !
1478 ! get the mean cooling in the range of z=z1 to z=z2
1479 !
1480  dtc = zero
1481  if ( zc > zero ) then
1482  if ( z1 < z2) then
1483  if ( z2 < zc ) then
1484  dtc = dt_cool*(one-(z1+z2)/(zc+zc))
1485  elseif ( z1 < zc .and. z2 >= zc ) then
1486  dtc = half*(one-z1/zc)*dt_cool*(zc-z1)/(z2-z1)
1487  endif
1488  elseif ( z1 == z2 ) then
1489  if ( z1 < zc ) then
1490  dtc = dt_cool*(one-z1/zc)
1491  endif
1492  endif
1493  endif
1494 
1495 !
1496 ! get the mean T departure from Tf in the range of z=z1 to z=z2
1497 !
1498  dtzm = dtw - dtc
1499 
1500  END SUBROUTINE dtzm_point
1501 
1521  subroutine tf_thaw_set(tf_ij,mask_ij,itile,jtile,tice,tclm,tf_thaw,nx,ny, &
1522  nset_thaw_s,nset_thaw_i,nset_thaw_c)
1523 
1524  real, dimension(nx*ny), intent(in) :: tf_ij
1525  integer, dimension(nx*ny), intent(in) :: mask_ij
1526  real, intent(in) :: tice,tclm
1527  integer, intent(in) :: itile,jtile,nx,ny
1528  real, intent(out) :: tf_thaw
1529  integer, intent(inout) :: nset_thaw_s,nset_thaw_i,nset_thaw_c
1530 ! Local
1531  real, parameter :: bmiss = -999.0
1532  real, dimension(nx,ny) :: tf
1533  integer, dimension(nx,ny) :: mask
1534  integer :: krad,max_search,istart,iend,jstart,jend
1535  integer :: ii,jj,iii,jjj
1536  logical :: is_ice
1537 
1538  max_search = 2
1539 
1540  mask(:,:) = reshape(mask_ij,(/nx,ny/) )
1541  tf(:,:) = reshape(tf_ij,(/nx,ny/) )
1542 
1543  tf_thaw = bmiss
1544 
1545  do krad = 1, max_search
1546 
1547  istart = itile - krad
1548  iend = itile + krad
1549  jstart = jtile - krad
1550  jend = jtile + krad
1551 
1552  do jj = jstart, jend
1553  do ii = istart, iend
1554 
1555 
1556  if ((jj == jstart) .or. (jj == jend) .or. &
1557  (ii == istart) .or. (ii == iend)) then
1558 
1559  if ((jj >= 1) .and. (jj <= ny)) then
1560  jjj = jj
1561  if (ii <= 0) then
1562  iii = nx + ii
1563  else if (ii >= (nx+1)) then
1564  iii = ii - nx
1565  else
1566  iii = ii
1567  endif
1568 
1569 !----------------------------------------------------------------------
1570 ! SEE IF NEARBY POINTS ARE SEA ICE. IF THEY ARE, AND THE SEARCH FOR
1571 ! A GAUSSIAN GRID OPEN WATER POINT FAILS, THEN TREF WILL BE SET TO
1572 ! FREEZING BELOW.
1573 !----------------------------------------------------------------------
1574  if (krad <= 2 .and. mask(iii,jjj) == 2) is_ice = .true.
1575 
1576  if (mask(iii,jjj) == 0) then
1577  tf_thaw = tf(iii,jjj)
1578  nset_thaw_s = nset_thaw_s + 1
1579  write(*,'(a,I4,2F9.3)') 'nset_thaw_s,tf(iii,jjj),tclm : ',nset_thaw_s,tf(iii,jjj),tclm
1580  go to 100
1581  endif ! tile mask is open water
1582 
1583  endif
1584  endif
1585  enddo
1586  enddo
1587  enddo ! krad loop
1588 
1589  100 continue
1590 
1591  if ( tf_thaw == bmiss ) then
1592  if (is_ice) then
1593  tf_thaw = tice
1594  nset_thaw_i = nset_thaw_i + 1
1595  write(*,'(a,I4,F9.3)') 'nset_thaw_i,tf_ice : ',nset_thaw_i,tice
1596  else
1597  tf_thaw = 0.8*tice+0.2*tclm
1598  nset_thaw_c = nset_thaw_c + 1
1599  write(*,'(a,I4,2F9.3)') 'nset_thaw_c,tf_ice,tclm : ',nset_thaw_c,tice,tclm
1600  endif
1601  endif
1602 
1603  end subroutine tf_thaw_set
1604 
1612  subroutine nsst_water_reset(nsst,ij,tf_thaw)
1613  use read_write_data, only : nsst_data
1614  implicit none
1615 
1616  integer, intent(in) :: ij
1617 
1618  real, intent(in) :: tf_thaw
1619 
1620  type(nsst_data), intent(inout) :: nsst
1621 
1622  nsst%c_0(ij) = 0.0
1623  nsst%c_d(ij) = 0.0
1624  nsst%d_conv(ij) = 0.0
1625  nsst%dt_cool(ij) = 0.0
1626  nsst%ifd(ij) = 0.0
1627  nsst%qrain(ij) = 0.0
1628  nsst%tref(ij) = tf_thaw
1629  nsst%w_0(ij) = 0.0
1630  nsst%w_d(ij) = 0.0
1631  nsst%xs(ij) = 0.0
1632  nsst%xt(ij) = 0.0
1633  nsst%xtts(ij) = 0.0
1634  nsst%xu(ij) = 0.0
1635  nsst%xv(ij) = 0.0
1636  nsst%xz(ij) = 30.0
1637  nsst%xzts(ij) = 0.0
1638  nsst%z_c(ij) = 0.0
1639  nsst%zm(ij) = 0.0
1640 
1641  end subroutine nsst_water_reset
1642 
1658 subroutine get_tf_clm(xlats_ij,xlons_ij,ny,nx,iy,im,id,ih,tf_clm,tf_trd)
1659  use read_write_data, only : get_tf_clm_dim
1660 
1661  implicit none
1662 
1663  real, dimension(nx*ny), intent(in) :: xlats_ij
1664  real, dimension(nx*ny), intent(in) :: xlons_ij
1665  real, dimension(nx,ny), intent(out) :: tf_clm
1666  real, dimension(nx,ny), intent(out) :: tf_trd
1667  integer, intent(in) :: iy,im,id,ih,nx,ny
1668 ! local declare
1669  real, allocatable, dimension(:,:) :: tf_clm0 ! sst climatology at the valid time (nxc,nyc)
1670  real, allocatable, dimension(:,:) :: tf_trd0 ! 6-hourly sst climatology tendency valid at atime (nxc,nyc)
1671  real, allocatable, dimension(:) :: cxlats ! latitudes of sst climatology
1672  real, allocatable, dimension(:) :: cxlons ! longitudes of sst climatology
1673 
1674  real, dimension(nx*ny) :: tf_clm_ij ! sst climatology at target grids (nx*ny)
1675  real, dimension(nx*ny) :: tf_trd_ij ! 6-hourly sst climatology tendency
1676  real :: wei1,wei2
1677  integer :: nxc,nyc,mon1,mon2
1678  character (len=6), parameter :: fin_tf_clm='sstclm' ! sst climatology file name
1679 !
1680 ! get which two months used and their weights from atime
1681 !
1682  call get_tim_wei(iy,im,id,ih,mon1,mon2,wei1,wei2)
1683 !
1684 ! get the dimensions of the sst climatology & allocate the related arrays
1685 !
1686  call get_tf_clm_dim(fin_tf_clm,nyc,nxc)
1687  allocate( tf_clm0(nxc,nyc),tf_trd0(nxc,nyc),cxlats(nyc),cxlons(nxc) )
1688 !
1689 ! get tf_clm at the analysis time from monthly climatology & cxlats, cxlons
1690 !
1691  call get_tf_clm_ta(tf_clm0,tf_trd0,cxlats,cxlons,nyc,nxc,mon1,mon2,wei1,wei2)
1692 !
1693 ! get tf_clm (nx by ny lat/lon) valid at atime
1694 !
1695  if ( nx == nxc .and. ny == nyc ) then
1696  tf_clm(:,:) = tf_clm0(:,:)
1697  tf_trd(:,:) = tf_trd0(:,:)
1698 ! write(*,'(a,2f9.3)') 'same dimensions, tf_clm, min : ',minval(tf_clm),maxval(tf_clm)
1699  else
1700 ! write(*,'(a,4i8)') 'different dimensions,nx,ny,nxc,nyc : ',nx,ny,nxc,nyc
1701  call intp_tile(tf_clm0, cxlats, cxlons, nyc, nxc, &
1702  tf_clm_ij,xlats_ij,xlons_ij,ny, nx)
1703  call intp_tile(tf_trd0, cxlats, cxlons, nyc, nxc, &
1704  tf_trd_ij,xlats_ij,xlons_ij,ny, nx)
1705 ! write(*,'(a,2f9.3)') 'tf_clm0, min, max : ',minval(tf_clm0),maxval(tf_clm0)
1706 
1707  tf_clm(:,:) = reshape(tf_clm_ij, (/nx,ny/) )
1708  tf_trd(:,:) = reshape(tf_trd_ij, (/nx,ny/) )
1709  endif
1710 
1711 end subroutine get_tf_clm
1712 
1727 subroutine get_tf_clm_ta(tf_clm_ta,tf_clm_trend,xlats,xlons,nlat,nlon,mon1,mon2,wei1,wei2)
1728  use read_write_data, only : read_tf_clim_grb
1729  implicit none
1730 
1731 ! input
1732  integer, intent(in) :: nlat,nlon,mon1,mon2
1733  real, intent(in) :: wei1,wei2
1734 ! output
1735  real, dimension(nlon,nlat), intent(out) :: tf_clm_ta,tf_clm_trend
1736  real, dimension(nlat), intent(out) :: xlats
1737  real, dimension(nlon), intent(out) :: xlons
1738 
1739 !input/output data file names
1740  character (len=6), parameter :: fin_tf_clm='sstclm'
1741 
1742 ! local declare
1743  real, dimension(nlon,nlat) :: tf_clm1,tf_clm2
1744 
1745 !
1746 ! read in rtg sst climatology without bitmap (surface mask) for mon1 and mon2
1747 !
1748  call read_tf_clim_grb(trim(fin_tf_clm),tf_clm1,xlats,xlons,nlat,nlon,mon1)
1749  call read_tf_clim_grb(trim(fin_tf_clm),tf_clm2,xlats,xlons,nlat,nlon,mon2)
1750 !
1751 ! tf_clim at the analysis time
1752 !
1753  tf_clm_ta(:,:) = wei1*tf_clm1(:,:)+wei2*tf_clm2(:,:)
1754 !
1755 ! tf_clim trend at the analysis time (month to month)
1756 !
1757  tf_clm_trend(:,:) = (tf_clm2(:,:)-tf_clm1(:,:))/120.0
1758 
1759  write(*,'(a,2f9.3)') 'tf_clm_ta, min, max : ',minval(tf_clm_ta),maxval(tf_clm_ta)
1760  write(*,'(a,2f9.3)') 'tf_clm_trend, min, max : ',minval(tf_clm_trend),maxval(tf_clm_trend)
1761  end subroutine get_tf_clm_ta
1762 
1775 subroutine get_sal_clm(xlats_ij,xlons_ij,ny,nx,iy,im,id,ih,sal_clm)
1776  use read_write_data, only : get_dim_nc
1777  implicit none
1778 
1779  real, dimension(nx*ny), intent(in) :: xlats_ij !
1780  real, dimension(nx*ny), intent(in) :: xlons_ij !
1781  real, dimension(nx,ny), intent(out) :: sal_clm !
1782  integer, intent(in) :: iy,im,id,ih,nx,ny
1783 ! local declare
1784  real, allocatable, dimension(:,:) :: sal_clm0 ! salinity climatology at the valid time
1785  real, allocatable, dimension(:) :: cxlats ! latitudes of sst climatology
1786  real, allocatable, dimension(:) :: cxlons ! longitudes of sst climatology
1787 
1788  real, dimension(nx*ny) :: sal_clm_ij ! salinity climatology at target grids (nx*ny)
1789  real :: wei1,wei2
1790  integer :: nxc,nyc,mon1,mon2
1791  character (len=6), parameter :: fin_sal_clm='salclm' ! salinity climatology file name
1792 !
1793 ! get which two months used and their weights from atime
1794 !
1795  call get_tim_wei(iy,im,id,ih,mon1,mon2,wei1,wei2)
1796 !
1797 ! get the dimensions of the sst climatology & allocate the related arrays
1798 !
1799  call get_dim_nc(fin_sal_clm,nyc,nxc)
1800  allocate( sal_clm0(nxc,nyc),cxlats(nyc),cxlons(nxc) )
1801 !
1802 ! get sal_clm at the analysis time from monthly climatology & cxlats, cxlons
1803 !
1804  call get_sal_clm_ta(sal_clm0,cxlats,cxlons,nyc,nxc,mon1,mon2,wei1,wei2)
1805 !
1806 ! get sal_clm (nx by ny lat/lon) valid at atime
1807 !
1808  if ( nx == nxc .and. ny == nyc ) then
1809  sal_clm(:,:) = sal_clm0(:,:)
1810 ! write(*,'(a,2f9.3)') 'same dimensions, sal_clm, min : ',minval(sal_clm),maxval(sal_clm)
1811  else
1812 ! write(*,'(a,4i8)') 'different dimensions,nx,ny,nxc,nyc : ',nx,ny,nxc,nyc
1813  call intp_tile(sal_clm0, cxlats, cxlons, nyc, nxc, &
1814  sal_clm_ij,xlats_ij,xlons_ij,ny, nx)
1815 ! write(*,'(a,2f9.3)') 'sal_clm0, min, max : ',minval(sal_clm0),maxval(sal_clm0)
1816 ! write(*,'(a,2f9.3)') 'done with intp_tile for sal_clm, min, max : ',minval(sal_clm_ij),maxval(sal_clm_ij)
1817 
1818  sal_clm(:,:) = reshape(sal_clm_ij, (/nx,ny/) )
1819  endif
1820 
1821 end subroutine get_sal_clm
1822 
1835 subroutine get_sal_clm_ta(sal_clm_ta,xlats,xlons,nlat,nlon,mon1,mon2,wei1,wei2)
1836 
1837  use read_write_data, only : read_salclm_gfs_nc
1838  implicit none
1839 
1840 ! input
1841  integer, intent(in) :: nlat,nlon,mon1,mon2
1842  real, intent(in) :: wei1,wei2
1843 ! output
1844  real, dimension(nlon,nlat), intent(out) :: sal_clm_ta
1845  real, dimension(nlat), intent(out) :: xlats
1846  real, dimension(nlon), intent(out) :: xlons
1847 
1848 !input/output data file names
1849  character (len=6), parameter :: fin_sal_clm='salclm'
1850 
1851 ! local declare
1852  real, dimension(nlon,nlat) :: sal_clm1,sal_clm2
1853 
1854 !
1855 ! read in rtg sst climatology without bitmap (surface mask) for mon1 and mon2
1856 !
1857  call read_salclm_gfs_nc(trim(fin_sal_clm),sal_clm1,xlats,xlons,nlat,nlon,mon1)
1858  call read_salclm_gfs_nc(trim(fin_sal_clm),sal_clm2,xlats,xlons,nlat,nlon,mon2)
1859 !
1860 ! sal_clim at the analysis time
1861 !
1862  sal_clm_ta(:,:) = wei1*sal_clm1(:,:)+wei2*sal_clm2(:,:)
1863  write(*,'(a,2f9.3)') 'sal_clm_ta, min, max : ',minval(sal_clm_ta),maxval(sal_clm_ta)
1864  end subroutine get_sal_clm_ta
1865 
1880 subroutine intp_tile(tf_lalo,dlats_lalo,dlons_lalo,jdim_lalo,idim_lalo, &
1881  tf_tile,xlats_tile,xlons_tile,jdim_tile,idim_tile)
1882 
1883  use utils
1884 
1885  implicit none
1886 
1887 ! input/output
1888  real, dimension(idim_lalo,jdim_lalo), intent(in) :: tf_lalo
1889  real, dimension(jdim_lalo), intent(in) :: dlats_lalo
1890  real, dimension(idim_lalo), intent(in) :: dlons_lalo
1891  real, dimension(jdim_tile*idim_tile), intent(in) :: xlats_tile
1892  real, dimension(jdim_tile*idim_tile), intent(in) :: xlons_tile
1893  integer, intent(in) :: jdim_lalo,idim_lalo,jdim_tile,idim_tile
1894  real, dimension(jdim_tile*idim_tile), intent(out) :: tf_tile
1895 
1896 ! local
1897  real, parameter :: deg2rad=3.1415926/180.0
1898  real, dimension(jdim_lalo) :: xlats_lalo
1899  real, dimension(idim_lalo) :: xlons_lalo
1900  real :: wsum
1901  integer :: itile,jtile
1902  integer :: ij
1903  integer :: ilalo,jlalo,ilalop1,jlalop1
1904 
1905  integer, allocatable, dimension(:,:) :: id1,id2,jdc
1906  real, allocatable, dimension(:,:,:) :: agrid,s2c
1907 
1908  print*
1909  print*,'interpolate from lat/lon grids to any one grid with known lat/lon'
1910 
1911  xlats_lalo = dlats_lalo*deg2rad
1912  xlons_lalo = dlons_lalo*deg2rad
1913 
1914  allocate(agrid(idim_tile,jdim_tile,2))
1915  agrid(:,:,1) = reshape(xlons_tile, (/idim_tile,jdim_tile/) )
1916  agrid(:,:,2) = reshape(xlats_tile, (/idim_tile,jdim_tile/) )
1917  agrid = agrid*deg2rad
1918 
1919  allocate(id1(idim_tile,jdim_tile))
1920  allocate(id2(idim_tile,jdim_tile))
1921  allocate(jdc(idim_tile,jdim_tile))
1922  allocate(s2c(idim_tile,jdim_tile,4))
1923 
1924 !----------------------------------------------------------------------
1925 ! compute bilinear weights for each model point from the nearest
1926 ! four lalo points. does not account for mask. that
1927 ! happens later.
1928 !----------------------------------------------------------------------
1929 
1930  call remap_coef( 1, idim_tile, 1, jdim_tile, idim_lalo, jdim_lalo, &
1931  xlons_lalo, xlats_lalo, id1, id2, jdc, s2c, agrid )
1932 
1933  do ij = 1, jdim_tile*idim_tile
1934 
1935  jtile = (ij-1)/idim_tile + 1
1936  itile = mod(ij,idim_tile)
1937  if (itile==0) itile = idim_tile
1938 
1939  ilalo = id1(itile,jtile)
1940  ilalop1 = id2(itile,jtile)
1941  jlalo = jdc(itile,jtile)
1942  jlalop1 = jdc(itile,jtile) + 1
1943 
1944  wsum = s2c(itile,jtile,1) + s2c(itile,jtile,2) + &
1945  s2c(itile,jtile,3) + s2c(itile,jtile,4)
1946 
1947  tf_tile(ij) = ( s2c(itile,jtile,1)*tf_lalo(ilalo,jlalo) + &
1948  s2c(itile,jtile,2)*tf_lalo(ilalop1,jlalo) + &
1949  s2c(itile,jtile,3)*tf_lalo(ilalop1,jlalop1) + &
1950  s2c(itile,jtile,4)*tf_lalo(ilalo,jlalop1) )/wsum
1951  enddo
1952 
1953  deallocate(id1, id2, jdc, s2c)
1954 
1955 end subroutine intp_tile
1956 
1969 subroutine get_tim_wei(iy,im,id,ih,mon1,mon2,wei1,wei2)
1970  implicit none
1971 
1972 ! input
1973  integer, intent(in) :: iy,im,id,ih
1974 ! output
1975  integer, intent(out) :: mon1,mon2
1976  real, intent(out) :: wei1,wei2
1977 
1978 ! local declare
1979  real :: rjday
1980  integer :: mon,monend,monm,monp,jdow,jdoy,jday
1981  integer :: jda(8)
1982 !
1983 !dayhf : julian day of the middle of each month
1984 !
1985  real, dimension(13) :: dayhf
1986  data dayhf/15.5,45.0,74.5,105.0,135.5,166.0,196.5,227.5,258.0,288.5,319.0,349.5,380.5/
1987 
1988 ! 15, 44, 73.5, 104, 134.5, 165, 195.5, 226.5, 257, 287.5, 318.5, 349 ! from
1989 ! woa05
1990 
1991  jda=0
1992  jda(1)=iy
1993  jda(2)=im
1994  jda(3)=id
1995  jda(5)=ih
1996  jdow = 0
1997  jdoy = 0
1998  jday = 0
1999  call w3doxdat(jda,jdow,jdoy,jday)
2000  rjday=jdoy+jda(5)/24.
2001  if(rjday.lt.dayhf(1)) rjday=rjday+365.
2002 
2003  monend = 12
2004  do mon = 1, monend
2005  monm = mon
2006  monp = mon + 1
2007  if( rjday >= dayhf(monm) .and. rjday < dayhf(monp) ) then
2008  mon1 = monm
2009  mon2 = monp
2010  go to 10
2011  endif
2012  enddo
2013 
2014  print *,'FATAL ERROR in get_tim_wei, wrong rjday',rjday
2015  call abort
2016  10 continue
2017 
2018  wei1 = (dayhf(mon2)-rjday)/(dayhf(mon2)-dayhf(mon1))
2019  wei2 = (rjday-dayhf(mon1))/(dayhf(mon2)-dayhf(mon1))
2020 
2021  if( mon2 == 13 ) mon2=1
2022 
2023  write(*,'(a,2i4,3f9.3)') 'mon1,mon2,rjday,wei1,wei2=',mon1,mon2,rjday,wei1,wei2
2024 
2025  end subroutine get_tim_wei
2026 
2036  real function tfreez(salinity)
2037 
2038  implicit none
2039 
2040  REAL salinity,sal
2041  REAL a1, a2, a3
2042  parameter(a1 = -0.0575)
2043  parameter(a2 = 1.710523e-3)
2044  parameter(a3 = -2.154996e-4)
2045 
2046  IF (salinity .LT. 0.) THEN
2047  sal = 0.
2048  ELSE
2049  sal = salinity
2050  ENDIF
2051  tfreez = sal*(a1+a2*sqrt(sal)+a3*sal)
2052 
2053  return
2054  end
num_parthds
integer function num_parthds()
Return the number of omp threads.
Definition: num_parthds.f90:9
land_increments::add_increment_soil
subroutine, public add_increment_soil(rla, rlo, stc_state, smc_state, slc_state, stc_updated, slc_updated, soilsnow_tile, soilsnow_fg_tile, lensfc, lsoil, idim, jdim, lsm, myrank)
Read in gsi file with soil state increments (on the gaussian grid), interpolate increments to the cub...
Definition: land_increments.f90:55
get_tf_clm
subroutine get_tf_clm(xlats_ij, xlons_ij, ny, nx, iy, im, id, ih, tf_clm, tf_trd)
Get the sst climatology at the valid time and on the target grid.
Definition: cycle.f90:1658
read_write_data::read_lat_lon_orog
subroutine, public read_lat_lon_orog(RLA, RLO, OROG, OROG_UF, TILE_NUM, IDIM, JDIM, IJDIM, LANDFRAC)
Read latitude and longitude for the cubed-sphere tile from the &#39;grid&#39; file.
Definition: read_write_data.f90:680
get_tf_clm_ta
subroutine get_tf_clm_ta(tf_clm_ta, tf_clm_trend, xlats, xlons, nlat, nlon, mon1, mon2, wei1, wei2)
Get the reference temperature/sst climatology and its trend at analysis time.
Definition: cycle.f90:1727
dtzm_point
subroutine dtzm_point(XT, XZ, DT_COOL, ZC, Z1, Z2, DTZM)
Compute the vertical mean of the NSST t-profile.
Definition: cycle.f90:1449
get_sal_clm_ta
subroutine get_sal_clm_ta(sal_clm_ta, xlats, xlons, nlat, nlon, mon1, mon2, wei1, wei2)
Get climatological salinity at the analysis time.
Definition: cycle.f90:1835
sfcdrv
subroutine sfcdrv(LUGB, IDIM, JDIM, LENSFC, LSOIL, DELTSFC, IY, IM, ID, IH, FH, IALB, USE_UFO, DO_NSST, DO_SFCCYCLE, DO_LNDINC, FRAC_GRID, ZSEA1, ZSEA2, ISOT, IVEGSRC, MYRANK)
Driver routine for updating the surface fields.
Definition: cycle.f90:299
tf_thaw_set
subroutine tf_thaw_set(tf_ij, mask_ij, itile, jtile, tice, tclm, tf_thaw, nx, ny, nset_thaw_s, nset_thaw_i, nset_thaw_c)
Set the background reference temperature (tf) for points where the ice has just melted.
Definition: cycle.f90:1521
land_increments::add_increment_snow
subroutine, public add_increment_snow(snd_inc, mask, lensfc, snd)
Add snow depth increment to model snow depth state, and limit output to be non-negative.
Definition: land_increments.f90:378
land_increments::calculate_landinc_mask
subroutine, public calculate_landinc_mask(smc, swe, vtype, lensfc, veg_type_landice, mask)
Calculate soil mask for land on model grid.
Definition: land_increments.f90:408
land_increments::apply_land_da_adjustments_snd
subroutine, public apply_land_da_adjustments_snd(lsm, lensfc, mask, swe_bck, snd_bck, snd_anl, swe_adj)
Make adjustments to dependent variables after applying land increments.
Definition: land_increments.f90:626
read_write_data::get_dim_nc
subroutine, public get_dim_nc(filename, nlat, nlon)
Get the i/j dimensions of the data from a NetCDF file.
Definition: read_write_data.f90:1824
machine
Holds machine dependent constants for global_cycle.
Definition: machine.f90:7
read_write_data::read_salclm_gfs_nc
subroutine, public read_salclm_gfs_nc(filename, sal, xlats, xlons, nlat, nlon, itime)
Read the woa05 salinity monthly climatology file.
Definition: read_write_data.f90:1753
land_increments
Definition: land_increments.f90:5
read_write_data::nsst_data
Definition: read_write_data.f90:16
intp_tile
subroutine intp_tile(tf_lalo, dlats_lalo, dlons_lalo, jdim_lalo, idim_lalo, tf_tile, xlats_tile, xlons_tile, jdim_tile, idim_tile)
Interpolate lon/lat grid data to the fv3 native grid (tf_lalo =&gt; tf_tile).
Definition: cycle.f90:1880
land_increments::apply_land_da_adjustments_soil
subroutine, public apply_land_da_adjustments_soil(lsm, isot, ivegsrc, lensfc, lsoil, rsoiltype, mask, stc_bck, stc_adj, smc_adj, slc_adj, stc_updated, slc_updated, zsoil)
Make adjustments to dependent variables after applying land increments.
Definition: land_increments.f90:465
climo_trend
subroutine climo_trend(LATITUDE, MON, DAY, DELTSFC, DTREF)
If the tile point is an isolated water point that has no corresponding gsi water point, then tref is updated using the rtg sst climo trend.
Definition: cycle.f90:1296
read_write_data::read_data
subroutine, public read_data(LSOIL, LENSFC, DO_NSST, INC_FILE, IS_NOAHMP, TSFFCS, SMCFCS, SWEFCS, STCFCS, TG3FCS, ZORFCS, CVFCS, CVBFCS, CVTFCS, ALBFCS, VEGFCS, SLIFCS, CNPFCS, F10M, VETFCS, SOTFCS, ALFFCS, USTAR, FMM, FHH, SIHFCS, SICFCS, SITFCS, TPRCP, SRFLAG, SNDFCS, VMNFCS, VMXFCS, SLCFCS, SLPFCS, ABSFCS, T2M, Q2M, SLMASK, ZSOIL, NSST)
Read the first guess surface records and nsst records (if selected) for a single cubed-sphere tile...
Definition: read_write_data.f90:1029
utils::remap_coef
subroutine, public remap_coef(is, ie, js, je, im, jm, lon, lat, id1, id2, jdc, s2c, agrid)
Generate the weights and index of the grids used in the bilinear interpolation.
Definition: utils.F90:38
read_write_data::get_tf_clm_dim
subroutine, public get_tf_clm_dim(file_sst, mlat_sst, mlon_sst)
Get the i/j dimensions of RTG SST climatology file.
Definition: read_write_data.f90:1695
sfc_drv
program sfc_drv
Stand alone surface/NSST cycle driver for the cubed-sphere grid.
Definition: cycle.f90:101
adjust_nsst
subroutine adjust_nsst(RLA, RLO, SLMSK_TILE, SLMSK_FG_TILE, SKINT_TILE, SICET_TILE, sice_tile, sice_fg_tile, SOILT_TILE, NSST, LENSFC, LSOIL, IDIM, JDIM, ZSEA1, ZSEA2, tf_clm_tile, tf_trd_tile, sal_clm_tile, LANDFRAC, FRAC_GRID)
Read in gsi file with the updated reference temperature increments (on the gaussian grid)...
Definition: cycle.f90:868
read_write_data::read_gsi_data
subroutine, public read_gsi_data(GSI_FILE, FILE_TYPE, LSOIL)
Read increment file from the GSI containing either the foundation temperature increments and mask...
Definition: read_write_data.f90:857
tfreez
real function tfreez(salinity)
Compute the freezing point of water as a function of salinity.
Definition: cycle.f90:2036
read_write_data
This module contains routines that read and write data.
Definition: read_write_data.f90:8
get_tim_wei
subroutine get_tim_wei(iy, im, id, ih, mon1, mon2, wei1, wei2)
For a given date, determine the bounding months and the linear time interpolation weights...
Definition: cycle.f90:1969
get_sal_clm
subroutine get_sal_clm(xlats_ij, xlons_ij, ny, nx, iy, im, id, ih, sal_clm)
Get salinity climatology at the valid time on the target grid.
Definition: cycle.f90:1775
nsst_water_reset
subroutine nsst_water_reset(nsst, ij, tf_thaw)
If the first guess was sea ice, but the analysis is open water, reset all nsst variables.
Definition: cycle.f90:1612
read_write_data::write_data
subroutine, public write_data(lensfc, idim, jdim, lsoil, do_nsst, nsst, slifcs, tsffcs, vegfcs, swefcs, tg3fcs, zorfcs, albfcs, alffcs, cnpfcs, f10m, t2m, q2m, vetfcs, sotfcs, ustar, fmm, fhh, sicfcs, sihfcs, sitfcs, tprcp, srflag, swdfcs, vmnfcs, vmxfcs, slpfcs, absfcs, slcfcs, smcfcs, stcfcs)
Update surface records - and nsst records if selected - on a single cubed-sphere tile to a pre-existi...
Definition: read_write_data.f90:120
utils
Module containing utility routines.
Definition: utils.F90:7
read_write_data::read_tf_clim_grb
subroutine, public read_tf_clim_grb(file_sst, sst, rlats_sst, rlons_sst, mlat_sst, mlon_sst, mon)
Read a GRIB1 sst climatological analysis file.
Definition: read_write_data.f90:1553