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