Changeset 10181


Ignore:
Timestamp:
2018-10-09T11:29:47+02:00 (20 months ago)
Author:
emmafiedler
Message:

Working version of ice thickness assimilation updates

Location:
branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/CONFIG/SHARED/namelist_ref

    r9987 r10181  
    12281228   ln_sla     = .false.             ! Logical switch for SLA observations 
    12291229   ln_sst     = .false.             ! Logical switch for SST observations 
    1230    ln_sic     = .false.             ! Logical switch for Sea Ice observations 
     1230   ln_sic     = .false.             ! Logical switch for Sea Ice concentration observations 
     1231   ln_sit     = .false.             ! Logical switch for Sea Ice thickness observations    
    12311232   ln_vel3d   = .false.             ! Logical switch for velocity observations 
    12321233   ln_sss     = .false.             ! Logical swithc for SSS observations 
     
    12671268   ln_sss_fp_indegs = .true. 
    12681269   ln_sic_fp_indegs = .true. 
     1270   ln_sit_fp_indegs = .true.    
    12691271! All of the *files* variables below are arrays. Use namelist_cfg to add more files 
    12701272   cn_profbfiles = 'profiles_01.nc'      ! Profile feedback input observation file names 
     
    12721274   cn_sstfbfiles = 'sst_01.nc'           ! SST feedback input observation file names 
    12731275   cn_sicfbfiles = 'sic_01.nc'           ! SIC feedback input observation file names 
     1276   cn_sitfbfiles = 'sit_01.nc'           ! SIT feedback input observation file names    
    12741277   cn_velfbfiles = 'vel_01.nc'           ! Velocity feedback input observation file names 
    12751278   cn_sssfbfiles = 'sss_01.nc'           ! SSS feedback input observation file names 
     
    13111314   rn_sic_avglamscl = 0.                 ! E/W diameter of SIC observation footprint (metres/degrees) 
    13121315   rn_sic_avgphiscl = 0.                 ! N/S diameter of SIC observation footprint (metres/degrees) 
     1316   rn_sit_avglamscl = 0.                 ! E/W diameter of SIT observation footprint (metres/degrees) 
     1317   rn_sit_avgphiscl = 0.                 ! N/S diameter of SIT observation footprint (metres/degrees) 
    13131318   nn_1dint = 0                          ! Type of vertical interpolation method 
    13141319   nn_2dint_default = 0                  ! Default horizontal interpolation method 
     
    13171322   nn_2dint_sss = -1                     ! Horizontal interpolation method for SSS 
    13181323   nn_2dint_sic = -1                     ! Horizontal interpolation method for SIC 
     1324   nn_2dint_sit = -1                     ! Horizontal interpolation method for SIT    
    13191325   nn_msshc = 0                          ! MSSH correction scheme 
    13201326   rn_mdtcorr = 1.61                     ! MDT  correction 
     
    13321338    ln_asmdin = .false.    !  Logical switch for Direct Initialization (DI) 
    13331339    ln_asmiau = .false.    !  Logical switch for Incremental Analysis Updating (IAU) 
    1334     ln_seaiceinc = .false. !  Logical switch for applying sea ice increments 
     1340    ln_seaiceinc = .false. !  Logical switch for applying sea ice concentration increments *UPDATE TO SIC* 
     1341    ln_sitinc = .false.    !  Logical switch for applying sea ice thickness increments 
    13351342    ln_temnofreeze = .false. !  Logical to not add increments if temperature would fall below freezing 
    13361343    nitbkg    = 0          !  Timestep of background in [0,nitend-nit000-1] 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r9987 r10181  
    2121   !!   dyn_asm_inc    : Apply the dynamic (u and v) increments 
    2222   !!   ssh_asm_inc    : Apply the SSH increment 
    23    !!   seaice_asm_inc : Apply the seaice increment 
     23   !!   seaice_asm_inc : Apply the seaice concentration increment 
     24   !!   sit_asm_inc    : Apply the sea ice thickness increment 
    2425   !!---------------------------------------------------------------------- 
    2526   USE wrk_nemo         ! Memory Allocation 
     
    4142#if defined key_cice && defined key_asminc 
    4243   USE sbc_ice, ONLY : & ! CICE Ice model variables 
    43    & ndaice_da, nfresh_da, nfsalt_da 
     44   & ndaice_da, ndsit_da, nfresh_da, nfsalt_da 
    4445#endif 
    4546   USE sbc_oce          ! Surface boundary condition variables. 
     
    5354   PUBLIC   dyn_asm_inc    !: Apply the dynamic (u and v) increments 
    5455   PUBLIC   ssh_asm_inc    !: Apply the SSH increment 
    55    PUBLIC   seaice_asm_inc !: Apply the seaice increment 
     56   PUBLIC   seaice_asm_inc !: Apply the seaice concentration increment 
     57   PUBLIC   sit_asm_inc    !: Apply the seaice thickness increment 
    5658 
    5759#if defined key_asminc 
     
    6668   LOGICAL, PUBLIC :: ln_dyninc = .FALSE.      !: No dynamics (u and v) assimilation increments 
    6769   LOGICAL, PUBLIC :: ln_sshinc = .FALSE.      !: No sea surface height assimilation increment 
    68    LOGICAL, PUBLIC :: ln_seaiceinc             !: No sea ice concentration increment 
     70   LOGICAL, PUBLIC :: ln_seaiceinc = .FALSE.   !: No sea ice concentration increment 
     71   LOGICAL, PUBLIC :: ln_sitinc = .FALSE.      !: No sea ice thickness increment 
    6972   LOGICAL, PUBLIC :: ln_salfix = .FALSE.      !: Apply minimum salinity check 
    7073   LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing 
     
    8992   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment 
    9093   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc 
     94   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   sit_bkginc            ! Increment to the background sea ice thickness 
    9195 
    9296   !! * Substitutions 
     
    136140         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   & 
    137141         &                 ln_salfix, salfixmin, nn_divdmp,                & 
    138          &                 ln_seaiceinc, ln_temnofreeze 
     142         &                 ln_seaiceinc, ln_sitinc, ln_temnofreeze 
    139143      !!---------------------------------------------------------------------- 
    140144 
     
    142146      ! Read Namelist nam_asminc : assimilation increment interface 
    143147      !----------------------------------------------------------------------- 
     148      ln_sitinc = .FALSE. 
    144149      ln_seaiceinc = .FALSE. 
    145150      ln_temnofreeze = .FALSE. 
     
    165170         WRITE(numout,*) '      Logical switch for applying SSH increments               ln_sshinc = ', ln_sshinc 
    166171         WRITE(numout,*) '      Logical switch for Direct Initialization (DI)            ln_asmdin = ', ln_asmdin 
    167          WRITE(numout,*) '      Logical switch for applying sea ice increments        ln_seaiceinc = ', ln_seaiceinc 
     172         WRITE(numout,*) '      Logical switch for applying SIC increments               ln_seaiceinc = ', ln_seaiceinc 
     173         WRITE(numout,*) '      Logical switch for applying SIT increments               ln_sitinc = ', ln_sitinc 
    168174         WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau 
    169175         WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg 
     
    221227 
    222228      IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) & 
    223            .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) & 
    224          & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', & 
     229         &  .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) & 
     230         &  .OR.( ln_sitinc ).OR.( ln_seaiceinc ) )) & 
     231         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc,', & 
     232         &                ' ln_sitinc and ln_seaiceinc is set to .true.', & 
    225233         &                ' but ln_asmdin and ln_asmiau are both set to .false. :', & 
    226234         &                ' Inconsistent options') 
     
    230238         &                ' Type IAU weighting function is invalid') 
    231239 
    232       IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) & 
     240      IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ) &  
     241         & .AND.( .NOT. ln_sitinc ).AND.( .NOT. ln_seaiceinc ) & 
    233242         &                     )  & 
    234          & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', & 
     243         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_sitinc', & 
     244         &                ' and ln_seaiceinc are set to .false. :', & 
    235245         &                ' The assimilation increments are not applied') 
    236246 
     
    336346      ALLOCATE( ssh_bkginc(jpi,jpj)   ) 
    337347      ALLOCATE( seaice_bkginc(jpi,jpj)) 
     348      ALLOCATE( sit_bkginc(jpi,jpj)   ) 
    338349#if defined key_asminc 
    339350      ALLOCATE( ssh_iau(jpi,jpj)      ) 
     
    345356      ssh_bkginc(:,:) = 0.0 
    346357      seaice_bkginc(:,:) = 0.0 
     358      sit_bkginc(:,:) = 0.0       
    347359#if defined key_asminc 
    348360      ssh_iau(:,:)    = 0.0 
    349361#endif 
    350       IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN 
     362      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_sitinc ).OR.( ln_seaiceinc ) ) THEN 
    351363 
    352364         !-------------------------------------------------------------------- 
     
    411423            ! to allow for differences in masks 
    412424            WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0 
     425         ENDIF 
     426 
     427         IF ( ln_sitinc ) THEN 
     428            CALL iom_get( inum, jpdom_autoglo, 'bckinsit', sit_bkginc, 1 ) 
     429            ! Apply the masks 
     430            sit_bkginc(:,:) = sit_bkginc(:,:) * tmask(:,:,1) 
     431            ! Set missing increments to 0.0 rather than 1e+20 
     432            ! to allow for differences in masks 
     433            WHERE( ABS( sit_bkginc(:,:) ) > 1.0e+10 ) sit_bkginc(:,:) = 0.0 
    413434         ENDIF 
    414435 
     
    769790         !   
    770791      ENDIF 
    771       ! Perhaps the following call should be in step 
    772       IF   ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment 
     792      ! Perhaps the following call should be in step       
     793      IF   ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment   
     794      IF   ( ln_sitinc  )      CALL sit_asm_inc ( kt )   ! apply sea ice thickness increment 
    773795      ! 
    774796   END SUBROUTINE tra_asm_inc 
     
    932954   END SUBROUTINE ssh_asm_inc 
    933955 
    934  
     956   SUBROUTINE sit_asm_inc( kt, kindic ) 
     957      !!---------------------------------------------------------------------- 
     958      !!                    ***  ROUTINE sit_asm_inc  *** 
     959      !!           
     960      !! ** Purpose : Apply the sea ice thickness assimilation increment. 
     961      !! 
     962      !! ** Method  : Direct initialization or Incremental Analysis Updating. 
     963      !! 
     964      !! ** Action  :  
     965      !! 
     966      !!---------------------------------------------------------------------- 
     967      IMPLICIT NONE 
     968      ! 
     969      INTEGER, INTENT(in)           ::   kt   ! Current time step 
     970      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation 
     971      ! 
     972      INTEGER  ::   it 
     973      REAL(wp) ::   zincwgt   ! IAU weight for current time step 
     974! #if defined key_lim2 
     975!       REAL(wp), DIMENSION(jpi,jpj) ::   zofrld, zohicif, zseaicendg, zhicifinc  ! LIM 
     976!       REAL(wp) ::   zhicifmin = 0.5_wp      ! ice minimum depth in metres 
     977!       !!THICKNESS INCS NOT SET UP FOR LIM 
     978! #endif 
     979      !!---------------------------------------------------------------------- 
     980 
     981      IF ( ln_asmiau ) THEN 
     982 
     983         !-------------------------------------------------------------------- 
     984         ! Incremental Analysis Updating 
     985         !-------------------------------------------------------------------- 
     986 
     987         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 
     988 
     989            it = kt - nit000 + 1 
     990            zincwgt = wgtiau(it)      ! IAU weight for the current time step  
     991            ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments) 
     992 
     993            IF(lwp) THEN 
     994               WRITE(numout,*)  
     995               WRITE(numout,*) 'sit_asm_inc : sea ice thick IAU at time step = ', & 
     996                  &  kt,' with IAU weight = ', wgtiau(it) 
     997               WRITE(numout,*) '~~~~~~~~~~~~' 
     998            ENDIF 
     999 
     1000            ! Sea-ice : LIM-3 case (to add) 
     1001 
     1002! #if defined key_lim2 
     1003!             ! Sea-ice : LIM-2 case (to add if needed) 
     1004!             zofrld (:,:) = frld(:,:) 
     1005!             zohicif(:,:) = hicif(:,:) 
     1006!             ! 
     1007!             frld  = MIN( MAX( frld (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp) 
     1008!             pfrld = MIN( MAX( pfrld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp) 
     1009!             fr_i(:,:) = 1.0_wp - frld(:,:)        ! adjust ice fraction 
     1010!             ! 
     1011!             zseaicendg(:,:) = zofrld(:,:) - frld(:,:)   ! find out actual sea ice nudge applied 
     1012!             ! 
     1013!             ! Nudge sea ice depth to bring it up to a required minimum depth 
     1014!             WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin )  
     1015!                zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt     
     1016!             ELSEWHERE 
     1017!                zhicifinc(:,:) = 0.0_wp 
     1018!             END WHERE 
     1019!             ! 
     1020!             ! nudge ice depth 
     1021!             hicif (:,:) = hicif (:,:) + zhicifinc(:,:) 
     1022!             phicif(:,:) = phicif(:,:) + zhicifinc(:,:)        
     1023!             ! 
     1024!             ! seaice salinity balancing (to add) 
     1025! #endif 
     1026 
     1027#if defined key_cice && defined key_asminc 
     1028            ! Sea-ice thickness : CICE case. Pass ice thickness increment tendency into CICE 
     1029            ndsit_da(:,:) = sit_bkginc(:,:) * zincwgt / rdt 
     1030#endif 
     1031 
     1032            IF ( kt == nitiaufin_r ) THEN 
     1033               DEALLOCATE( sit_bkginc ) 
     1034            ENDIF 
     1035 
     1036         ELSE 
     1037 
     1038#if defined key_cice && defined key_asminc 
     1039            ! Sea-ice thickness : CICE case. Zero ice increment tendency into CICE 
     1040            ndsit_da(:,:) = 0.0_wp 
     1041#endif 
     1042 
     1043         ENDIF 
     1044 
     1045      ELSEIF ( ln_asmdin ) THEN 
     1046 
     1047         !-------------------------------------------------------------------- 
     1048         ! Direct Initialization 
     1049         !-------------------------------------------------------------------- 
     1050 
     1051         IF ( kt == nitdin_r ) THEN 
     1052 
     1053            neuler = 0                    ! Force Euler forward step 
     1054 
     1055            ! Sea-ice : LIM-3 case (to add) 
     1056 
     1057! #if defined key_lim2 
     1058!             ! Sea-ice : LIM-2 case (add if needed) 
     1059!             zofrld(:,:)=frld(:,:) 
     1060!             zohicif(:,:)=hicif(:,:) 
     1061!             !  
     1062!             ! Initialize the now fields the background + increment 
     1063!             frld (:,:) = MIN( MAX( frld(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp) 
     1064!             pfrld(:,:) = frld(:,:)  
     1065!             fr_i (:,:) = 1.0_wp - frld(:,:)                ! adjust ice fraction 
     1066!             zseaicendg(:,:) = zofrld(:,:) - frld(:,:)      ! find out actual sea ice nudge applied 
     1067!             ! 
     1068!             ! Nudge sea ice depth to bring it up to a required minimum depth 
     1069!             WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin )  
     1070!                zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt     
     1071!             ELSEWHERE 
     1072!                zhicifinc(:,:) = 0.0_wp 
     1073!             END WHERE 
     1074!             ! 
     1075!             ! nudge ice depth 
     1076!             hicif (:,:) = hicif (:,:) + zhicifinc(:,:) 
     1077!             phicif(:,:) = phicif(:,:)        
     1078!             ! 
     1079!             ! seaice salinity balancing (to add) 
     1080! #endif 
     1081  
     1082#if defined key_cice && defined key_asminc 
     1083            ! Sea-ice thickness : CICE case. Pass ice thickness increment tendency into CICE 
     1084           ndsit_da(:,:) = sit_bkginc(:,:) / rdt 
     1085#endif 
     1086           IF ( .NOT. PRESENT(kindic) ) THEN 
     1087              DEALLOCATE( sit_bkginc ) 
     1088           END IF 
     1089 
     1090         ELSE 
     1091 
     1092#if defined key_cice && defined key_asminc 
     1093            ! Sea-ice thicnkness : CICE case. Zero ice thickness increment tendency into CICE  
     1094            ndsit_da(:,:) = 0.0_wp 
     1095#endif 
     1096          
     1097         ENDIF 
     1098 
     1099!#if defined defined key_lim2 || defined key_cice 
     1100! 
     1101!            IF (ln_seaicebal ) THEN        
     1102!             !! balancing salinity increments 
     1103!             !! simple case from limflx.F90 (doesn't include a mass flux) 
     1104!             !! assumption is that as ice concentration is reduced or increased 
     1105!             !! the snow and ice depths remain constant 
     1106!             !! note that snow is being created where ice concentration is being increased 
     1107!             !! - could be more sophisticated and 
     1108!             !! not do this (but would need to alter h_snow) 
     1109! 
     1110!             usave(:,:,:)=sb(:,:,:)   ! use array as a temporary store 
     1111! 
     1112!             DO jj = 1, jpj 
     1113!               DO ji = 1, jpi  
     1114!           ! calculate change in ice and snow mass per unit area 
     1115!           ! positive values imply adding salt to the ocean (results from ice formation) 
     1116!           ! fwf : ice formation and melting 
     1117! 
     1118!                 zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rdt 
     1119! 
     1120!           ! change salinity down to mixed layer depth 
     1121!                 mld=hmld_kara(ji,jj) 
     1122! 
     1123!           ! prevent small mld 
     1124!           ! less than 10m can cause salinity instability  
     1125!                 IF (mld < 10) mld=10 
     1126! 
     1127!           ! set to bottom of a level  
     1128!                 DO jk = jpk-1, 2, -1 
     1129!                   IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN  
     1130!                     mld=gdepw(ji,jj,jk+1) 
     1131!                     jkmax=jk 
     1132!                   ENDIF 
     1133!                 ENDDO 
     1134! 
     1135!            ! avoid applying salinity balancing in shallow water or on land 
     1136!            !  
     1137! 
     1138!            ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m) 
     1139! 
     1140!                 dsal_ocn=0.0_wp 
     1141!                 sal_thresh=5.0_wp        ! minimum salinity threshold for salinity balancing 
     1142! 
     1143!                 if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) & 
     1144!                              dsal_ocn = zfons / (rhop(ji,jj,1) * mld) 
     1145! 
     1146!           ! put increments in for levels in the mixed layer 
     1147!           ! but prevent salinity below a threshold value  
     1148! 
     1149!                   DO jk = 1, jkmax               
     1150! 
     1151!                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN  
     1152!                           sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn 
     1153!                           sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn 
     1154!                     ENDIF 
     1155! 
     1156!                   ENDDO 
     1157! 
     1158!      !            !  salt exchanges at the ice/ocean interface 
     1159!      !            zpmess         = zfons / rdt_ice    ! rdt_ice is ice timestep 
     1160!      ! 
     1161!      !! Adjust fsalt. A +ve fsalt means adding salt to ocean 
     1162!      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt   
     1163!      !!                
     1164!      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d)  
     1165!      !!                                                     ! E-P (kg m-2 s-2) 
     1166!      !            emp(ji,jj) = emp(ji,jj) + zpmess          ! E-P (kg m-2 s-2) 
     1167!               ENDDO !ji 
     1168!             ENDDO !jj! 
     1169! 
     1170!            ENDIF !ln_seaicebal 
     1171! 
     1172!#endif 
     1173 
     1174      ENDIF 
     1175 
     1176   END SUBROUTINE sit_asm_inc 
     1177    
    9351178   SUBROUTINE seaice_asm_inc( kt, kindic ) 
    9361179      !!---------------------------------------------------------------------- 
    9371180      !!                    ***  ROUTINE seaice_asm_inc  *** 
    9381181      !!           
    939       !! ** Purpose : Apply the sea ice assimilation increment. 
     1182      !! ** Purpose : Apply the sea ice concentration assimilation increment. 
    9401183      !! 
    9411184      !! ** Method  : Direct initialization or Incremental Analysis Updating. 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

    r9306 r10181  
    5252   LOGICAL :: ln_sst_fp_indegs     !: T=>     SST obs footprint size specified in degrees, F=> in metres 
    5353   LOGICAL :: ln_sss_fp_indegs     !: T=>     SSS obs footprint size specified in degrees, F=> in metres 
    54    LOGICAL :: ln_sic_fp_indegs     !: T=> sea-ice obs footprint size specified in degrees, F=> in metres 
     54   LOGICAL :: ln_sic_fp_indegs     !: T=> SIC obs footprint size specified in degrees, F=> in metres 
     55   LOGICAL :: ln_sit_fp_indegs     !: T=> SIT obs footprint size specified in degrees, F=> in metres 
    5556 
    5657   REAL(wp) :: rn_default_avglamscl !: Default E/W diameter of observation footprint 
     
    6263   REAL(wp) :: rn_sss_avglamscl     !: E/W diameter of SSS observation footprint 
    6364   REAL(wp) :: rn_sss_avgphiscl     !: N/S diameter of SSS observation footprint 
    64    REAL(wp) :: rn_sic_avglamscl     !: E/W diameter of sea-ice observation footprint 
    65    REAL(wp) :: rn_sic_avgphiscl     !: N/S diameter of sea-ice observation footprint 
    66  
     65   REAL(wp) :: rn_sic_avglamscl     !: E/W diameter of SIC observation footprint 
     66   REAL(wp) :: rn_sic_avgphiscl     !: N/S diameter of SIC observation footprint 
     67   REAL(wp) :: rn_sit_avglamscl     !: E/W diameter of SIT observation footprint 
     68   REAL(wp) :: rn_sit_avgphiscl     !: N/S diameter of SIT observation footprint 
     69    
    6770   INTEGER :: nn_1dint         !: Vertical interpolation method 
    6871   INTEGER :: nn_2dint_default !: Default horizontal interpolation method 
     
    7073   INTEGER :: nn_2dint_sst     !: SST horizontal interpolation method (-1 = default) 
    7174   INTEGER :: nn_2dint_sss     !: SSS horizontal interpolation method (-1 = default) 
    72    INTEGER :: nn_2dint_sic     !: Seaice horizontal interpolation method (-1 = default) 
     75   INTEGER :: nn_2dint_sic     !: SIC horizontal interpolation method (-1 = default) 
     76   INTEGER :: nn_2dint_sit     !: SIT horizontal interpolation method (-1 = default) 
    7377  
    7478   INTEGER, DIMENSION(imaxavtypes) :: & 
     
    151155         & cn_slafbfiles,      & ! Sea level anomaly input filenames 
    152156         & cn_sicfbfiles,      & ! Seaice concentration input filenames 
     157         & cn_sitfbfiles,      & ! Seaice thickness input filenames          
    153158         & cn_velfbfiles,      & ! Velocity profile input filenames 
    154159         & cn_sssfbfiles,      & ! Sea surface salinity input filenames 
     
    187192      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature 
    188193      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration 
     194      LOGICAL :: ln_sit          ! Logical switch for sea ice thickness      
    189195      LOGICAL :: ln_sss          ! Logical switch for sea surface salinity obs 
    190196      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs 
     
    242248 
    243249      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              & 
    244          &            ln_sst, ln_sic, ln_sss, ln_vel3d,               & 
     250         &            ln_sst, ln_sic, ln_sit, ln_sss, ln_vel3d,       & 
    245251         &            ln_slchltot, ln_slchldia, ln_slchlnon,          & 
    246252         &            ln_slchldin, ln_slchlmic, ln_slchlnan,          & 
     
    257263         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             & 
    258264         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             & 
     265         &            ln_sit_fp_indegs,                               & 
    259266         &            cn_profbfiles, cn_slafbfiles,                   & 
    260267         &            cn_sstfbfiles, cn_sicfbfiles,                   & 
     268         &            cn_sitfbfiles,                                  & 
    261269         &            cn_velfbfiles, cn_sssfbfiles,                   & 
    262270         &            cn_slchltotfbfiles, cn_slchldiafbfiles,         & 
     
    279287         &            rn_sss_avglamscl, rn_sss_avgphiscl,             & 
    280288         &            rn_sic_avglamscl, rn_sic_avgphiscl,             & 
     289         &            rn_sit_avglamscl, rn_sit_avgphiscl,             & 
    281290         &            nn_1dint, nn_2dint_default,                     & 
    282291         &            nn_2dint_sla, nn_2dint_sst,                     & 
    283292         &            nn_2dint_sss, nn_2dint_sic,                     & 
     293         &            nn_2dint_sit,                                   & 
    284294         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             & 
    285295         &            nn_profdavtypes 
     
    294304      cn_sstfbfiles(:)      = '' 
    295305      cn_sicfbfiles(:)      = '' 
     306      cn_sitfbfiles(:)      = ''       
    296307      cn_velfbfiles(:)      = '' 
    297308      cn_sssfbfiles(:)      = '' 
     
    355366         WRITE(numout,*) '             Logical switch for SLA observations                      ln_sla = ', ln_sla 
    356367         WRITE(numout,*) '             Logical switch for SST observations                      ln_sst = ', ln_sst 
    357          WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic 
     368         WRITE(numout,*) '             Logical switch for Sea Ice Concentration observations    ln_sic = ', ln_sic 
     369         WRITE(numout,*) '             Logical switch for Sea Ice Thickness observations        ln_sit = ', ln_sit          
    358370         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d 
    359371         WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss 
     
    393405         WRITE(numout,*) '             Type of horizontal interpolation method for SSS    nn_2dint_sss = ', nn_2dint_sss 
    394406         WRITE(numout,*) '             Type of horizontal interpolation method for SIC    nn_2dint_sic = ', nn_2dint_sic 
     407         WRITE(numout,*) '             Type of horizontal interpolation method for SIT    nn_2dint_sit = ', nn_2dint_sit          
    395408         WRITE(numout,*) '             Default E/W diameter of obs footprint      rn_default_avglamscl = ', rn_default_avglamscl 
    396409         WRITE(numout,*) '             Default N/S diameter of obs footprint      rn_default_avgphiscl = ', rn_default_avgphiscl 
     
    405418         WRITE(numout,*) '             SIC N/S diameter of obs footprint              rn_sic_avgphiscl = ', rn_sic_avgphiscl 
    406419         WRITE(numout,*) '             SIC obs footprint in deg [T] or m [F]          ln_sic_fp_indegs = ', ln_sic_fp_indegs 
     420         WRITE(numout,*) '             SIT E/W diameter of obs footprint              rn_sit_avglamscl = ', rn_sit_avglamscl 
     421         WRITE(numout,*) '             SIT N/S diameter of obs footprint              rn_sit_avgphiscl = ', rn_sit_avgphiscl 
     422         WRITE(numout,*) '             SIT obs footprint in deg [T] or m [F]          ln_sit_fp_indegs = ', ln_sit_fp_indegs          
    407423         WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea 
    408424         WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject 
     
    424440         &                  ln_pchltot,  ln_pno3,     ln_psi4,     ln_ppo4,     & 
    425441         &                  ln_pdic,     ln_palk,     ln_pph,      ln_po2 /) ) 
    426       nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss,                     & 
     442      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sit, ln_sss,             & 
    427443         &                  ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, & 
    428444         &                  ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot,  & 
     
    535551            clsurffiles(jtype,:) = cn_sicfbfiles 
    536552         ENDIF 
     553         IF (ln_sit) THEN 
     554            jtype = jtype + 1 
     555            cobstypessurf(jtype) = 'sit' 
     556            clsurffiles(jtype,:) = cn_sitfbfiles 
     557         ENDIF          
    537558         IF (ln_sss) THEN 
    538559            jtype = jtype + 1 
     
    645666               ltype_fp_indegs = ln_sic_fp_indegs 
    646667               ltype_night     = .FALSE. 
     668            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sit' ) THEN 
     669               IF ( nn_2dint_sit == -1 ) THEN 
     670                  n2dint_type  = nn_2dint_default 
     671               ELSE 
     672                  n2dint_type  = nn_2dint_sit 
     673               ENDIF 
     674               ztype_avglamscl = rn_sit_avglamscl 
     675               ztype_avgphiscl = rn_sit_avgphiscl 
     676               ltype_fp_indegs = ln_sit_fp_indegs 
     677               ltype_night     = .FALSE.                
    647678            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 
    648679               IF ( nn_2dint_sss == -1 ) THEN 
     
    877908#endif 
    878909#if defined key_cice 
    879       USE sbc_oce, ONLY : fr_i     ! ice fraction 
     910      USE sbc_oce, ONLY : &        ! CICE variables 
     911         & fr_i,          &        ! ice fraction 
     912         & thick_i                 ! ice thickness 
    880913#endif 
    881914#if defined key_hadocc 
     
    11901223                        &           'time-step but some obs are valid then.' ) 
    11911224                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 
    1192                         &           ' sea-ice obs will be missed' 
     1225                        &           ' sea-ice concentration obs will be missed' 
    11931226                  ENDIF 
    11941227                  surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + & 
     
    12011234                  zsurfvar(:,:) = 1._wp - frld(:,:) 
    12021235#else 
    1203                CALL ctl_stop( ' Trying to run sea-ice observation operator', & 
     1236               CALL ctl_stop( ' Trying to run sea-ice concentration observation operator', & 
     1237                  &           ' but no sea-ice model appears to have been defined' ) 
     1238#endif 
     1239               ENDIF 
     1240            CASE('sit') 
     1241               IF ( kstp == 0 ) THEN ! **Copied from SIC, check applies to SIT!** 
     1242                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 
     1243                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 
     1244                        &           'time-step but some obs are valid then.' ) 
     1245                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 
     1246                        &           ' sea-ice thickness obs will be missed' 
     1247                  ENDIF 
     1248                  surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + & 
     1249                     &                        surfdataqc(jtype)%nsstp(1) 
     1250                  CYCLE 
     1251               ELSE             
     1252#if defined key_cice 
     1253               ! Will need to convert freeboard to thickness here 
     1254                  zsurfvar(:,:) = thick_i(:,:) 
     1255#elif defined key_lim2 || defined key_lim3 
     1256               CALL ctl_stop( ' No sea-ice thickness observation operator defined for LIM model' ) 
     1257#else 
     1258               CALL ctl_stop( ' Trying to run sea-ice thickness observation operator', & 
    12041259                  &           ' but no sea-ice model appears to have been defined' ) 
    12051260#endif 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90

    r9308 r10181  
    493493         clgrid     = 'T' 
    494494 
     495      CASE('SIT') 
     496 
     497         clfiletype = 'sitfb' 
     498         cllongname = 'Sea ice thickness' 
     499         clunits    = 'm' 
     500         clgrid     = 'T' 
     501          
    495502      CASE('SSS') 
    496503 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

    r9987 r10181  
    105105   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tsfc_ice           !: sea-ice surface skin temperature (on categories) 
    106106   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   kn_ice             !: sea-ice surface layer thermal conductivity (on cats) 
    107  
     107   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   thick_iu           !: ice thickness at NEMO U point 
     108   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   thick_iv           !: ice thickness at NEMO V point 
    108109   ! variables used in the coupled interface 
    109110   INTEGER , PUBLIC, PARAMETER ::   jpl = ncat 
     
    116117#if defined key_asminc 
    117118   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ndaice_da          !: NEMO fresh water flux to ocean due to data assim 
     119   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ndsit_da           !: NEMO ice thickness change from data assim (/second? These descriptions are jumbled?) 
    118120   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nfresh_da          !: NEMO salt flux to ocean due to data assim 
    119121   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nfsalt_da          !: NEMO ice concentration change/second from data assim 
     
    171173                wndj_ice(jpi,jpj)     , nfrzmlt(jpi,jpj)      , ss_iou(jpi,jpj)       , & 
    172174                ss_iov(jpi,jpj)       , fr_iu(jpi,jpj)        , fr_iv(jpi,jpj)        , & 
     175                thick_iu(jpi,jpj)     , thick_iv(jpi,jpj)     ,                         & 
     176                ht_i(jpi,jpj,ncat)    , ht_s(jpi,jpj,ncat)    ,                         &                 
    173177                a_i(jpi,jpj,ncat)     , topmelt(jpi,jpj,ncat) , botmelt(jpi,jpj,ncat) , & 
    174178#if defined key_asminc 
    175179                ndaice_da(jpi,jpj)    , nfresh_da(jpi,jpj)    , nfsalt_da(jpi,jpj)    , & 
     180                ndsit_da(jpi,jpj)    ,                                                  & 
    176181#endif 
    177182                sstfrz(jpi,jpj)       , STAT= ierr(1) ) 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r9987 r10181  
    121121   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sprecip           !: solid precipitation                          [Kg/m2/s] 
    122122   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fr_i              !: ice fraction = 1 - lead fraction      (between 0 to 1) 
     123   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   thick_i           !: ice thickness [m] 
     124   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vol_i             !: ice volume [m3]    
    123125#if defined key_cpl_carbon_cycle 
    124126   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   atm_co2           !: atmospheric pCO2                             [ppm] 
     
    196198         ! 
    197199      ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     & 
     200         & thick_i(jpi,jpj) ,  vol_i(jpi,jpj) ,                           &       
    198201#if defined key_cpl_carbon_cycle 
    199202         &      atm_co2(jpi,jpj) ,                                        & 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90

    r9987 r10181  
    5757                flatn_f,fsurfn_f,fcondtopn_f,                    & 
    5858#ifdef key_asminc 
    59                 daice_da,fresh_da,fsalt_da,                    & 
     59                daice_da,dsit_da,fresh_da,fsalt_da,             & 
    6060#endif 
    6161                uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl,   & 
     
    250250      CALL lbc_lnk ( fr_iu , 'U', 1. ) 
    251251      CALL lbc_lnk ( fr_iv , 'V', 1. ) 
     252       
     253!!*****JUST ADDED BUT NEEDS CONFIRMING!!**** 
     254! Snow and ice thickness 
     255      CALL cice2nemo(vice,thick_i,'T', 1. )  !!Do this for snow depth too? set thick_s. Duplicated below for embedded ice? 
     256 
     257! vice is volume per unit area of grid cell = thickness 
     258       
     259      IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN ! Confirm what this if loop is for 
     260         DO jl = 1,ncat 
     261            CALL cice2nemo(vsnon(:,:,jl,:),ht_s(:,:,jl),'T', 1. ) 
     262            CALL cice2nemo(vicen(:,:,jl,:),ht_i(:,:,jl),'T', 1. ) 
     263         ENDDO 
     264      ENDIF 
     265       
     266! T point to U point 
     267! T point to V point 
     268      thick_iu(:,:)=0.0 
     269      thick_iv(:,:)=0.0 
     270      DO jj=1,jpjm1 
     271         DO ji=1,jpim1 
     272            thick_iu(ji,jj)=0.5*(thick_i(ji,jj)+thick_i(ji+1,jj))*umask(ji,jj,1) 
     273            thick_iv(ji,jj)=0.5*(thick_i(ji,jj)+thick_i(ji,jj+1))*vmask(ji,jj,1) 
     274         ENDDO 
     275      ENDDO 
     276 
     277      CALL lbc_lnk ( thick_iu , 'U', 1. ) 
     278      CALL lbc_lnk ( thick_iv , 'V', 1. ) 
     279 
     280!!**********************************************************      
     281       
    252282 
    253283      !                                      ! embedded sea ice 
     
    310340      nfresh_da(:,:) = 0.0    
    311341      nfsalt_da(:,:) = 0.0    
    312       ndaice_da(:,:) = 0.0          
     342      ndaice_da(:,:) = 0.0 
     343      ndsit_da(:,:) = 0.0 
    313344#endif 
    314345      ! 
     
    469500      ztmp(:,:)=ndaice_da(:,:)*tmask(:,:,1) 
    470501      Call nemo2cice(ztmp,daice_da,'T', 1. ) 
     502!Ice thickness change (from assimilation) 
     503      ztmp(:,:)=ndsit_da(:,:)*tmask(:,:,1) 
     504      Call nemo2cice(ztmp,dsit_da,'T', 1. ) 
    471505#endif  
    472506 
Note: See TracChangeset for help on using the changeset viewer.