New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 10177 – NEMO

Changeset 10177


Ignore:
Timestamp:
2018-10-08T11:16:10+02:00 (6 years ago)
Author:
csanchez
Message:

Added changes to syncronize with fcm:nemo.xm/branches/UKMO/AMM15_v3_6_STABLE_package@9537

Location:
branches/UKMO/CCRS_NUS_MC_STABLE/NEMOGCM/NEMO/OPA_SRC
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • branches/UKMO/CCRS_NUS_MC_STABLE/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90

    r8058 r10177  
    5050   USE ice 
    5151#endif 
     52   USE asminc, ONLY: ln_avgbkg 
    5253   IMPLICIT NONE 
    5354   PRIVATE 
    5455    
    5556   PUBLIC   asm_bkg_wri   !: Write out the background state 
     57 
     58  !! * variables for calculating time means 
     59   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_tavg  , sn_tavg   
     60   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   un_tavg  , vn_tavg 
     61   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   avt_tavg 
     62#if defined key_zdfgls || key_zdftke 
     63   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   en_tavg 
     64#endif 
     65   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_tavg 
     66   REAL(wp),SAVE :: numtimes_tavg     ! No of times to average over 
    5667 
    5768   !!---------------------------------------------------------------------- 
    … …  
    8192      INTEGER :: inum          ! File unit number 
    8293      REAL(wp) :: zdate        ! Date 
     94      INTEGER :: ierror 
    8395      !!----------------------------------------------------------------------- 
    8496 
    85       !                                !------------------------------------------- 
    86       IF( kt == nitbkg_r ) THEN        ! Write out background at time step nitbkg_r 
    87          !                             !-----------------------------------======== 
     97      ! If creating an averaged assim bkg, initialise on first timestep 
     98      IF ( ln_avgbkg .AND. kt == ( nn_it000 - 1) ) THEN 
     99         ! Allocate memory  
     100         ALLOCATE( tn_tavg(jpi,jpj,jpk), STAT=ierror ) 
     101         IF( ierror > 0 ) THEN 
     102            CALL ctl_stop( 'asm_wri_bkg: unable to allocate tn_tavg' )   ;   RETURN 
     103         ENDIF 
     104         tn_tavg(:,:,:)=0 
     105         ALLOCATE( sn_tavg(jpi,jpj,jpk), STAT=ierror ) 
     106         IF( ierror > 0 ) THEN 
     107            CALL ctl_stop( 'asm_wri_bkg: unable to allocate sn_tavg' )   ;   RETURN 
     108         ENDIF 
     109         sn_tavg(:,:,:)=0 
     110         ALLOCATE( un_tavg(jpi,jpj,jpk), STAT=ierror ) 
     111         IF( ierror > 0 ) THEN 
     112            CALL ctl_stop( 'asm_wri_bkg: unable to allocate un_tavg' )   ;   RETURN 
     113         ENDIF 
     114         un_tavg(:,:,:)=0 
     115         ALLOCATE( vn_tavg(jpi,jpj,jpk), STAT=ierror ) 
     116         IF( ierror > 0 ) THEN 
     117            CALL ctl_stop( 'asm_wri_bkg: unable to allocate vn_tavg' )   ;   RETURN 
     118         ENDIF 
     119         vn_tavg(:,:,:)=0 
     120         ALLOCATE( sshn_tavg(jpi,jpj), STAT=ierror ) 
     121         IF( ierror > 0 ) THEN 
     122            CALL ctl_stop( 'asm_wri_bkg: unable to allocate sshn_tavg' )   ;   RETURN 
     123         ENDIF 
     124         sshn_tavg(:,:)=0 
     125#if defined key_zdftke 
     126         ALLOCATE( en_tavg(jpi,jpj,jpk), STAT=ierror ) 
     127         IF( ierror > 0 ) THEN 
     128            CALL ctl_stop( 'asm_wri_bkg: unable to allocate en_tavg' )   ;   RETURN 
     129         ENDIF 
     130         en_tavg(:,:,:)=0 
     131#endif 
     132         ALLOCATE( avt_tavg(jpi,jpj,jpk), STAT=ierror ) 
     133         IF( ierror > 0 ) THEN 
     134            CALL ctl_stop( 'asm_wri_bkg: unable to allocate avt_tavg' )   ;   RETURN 
     135         ENDIF 
     136         avt_tavg(:,:,:)=0 
     137          
     138         numtimes_tavg = REAL ( nitavgbkg_r -  nn_it000 + 1 ) 
     139      ENDIF    
     140 
     141      ! If creating an averaged assim bkg, sum the contribution every timestep 
     142      IF ( ln_avgbkg ) THEN  
     143         IF (lwp) THEN 
     144              WRITE(numout,*) 'asm_wri_bkg : Summing assim bkg fields at timestep ',kt 
     145              WRITE(numout,*) '~~~~~~~~~~~~ ' 
     146         ENDIF 
     147 
     148         tn_tavg(:,:,:)        = tn_tavg(:,:,:) + tsn(:,:,:,jp_tem) / numtimes_tavg 
     149         sn_tavg(:,:,:)        = sn_tavg(:,:,:) + tsn(:,:,:,jp_sal) / numtimes_tavg 
     150         sshn_tavg(:,:)        = sshn_tavg(:,:) + sshn (:,:) / numtimes_tavg 
     151         un_tavg(:,:,:)        = un_tavg(:,:,:) + un(:,:,:) / numtimes_tavg 
     152         vn_tavg(:,:,:)        = vn_tavg(:,:,:) + vn(:,:,:) / numtimes_tavg 
     153         avt_tavg(:,:,:)        = avt_tavg(:,:,:) + avt(:,:,:) / numtimes_tavg 
     154#if defined key_zdftke 
     155         en_tavg(:,:,:)       = en_tavg(:,:,:) + en(:,:,:) / numtimes_tavg 
     156#endif 
     157      ENDIF 
     158      
     159 
     160      ! Write out background at time step nitbkg_r or nitavgbkg_r 
     161      IF ( ( .NOT. ln_avgbkg .AND. (kt == nitbkg_r) ) .OR. & 
     162      &          ( ln_avgbkg .AND. (kt == nitavgbkg_r) ) ) THEN 
    88163         ! 
    89164         WRITE(cl_asmbkg, FMT='(A,".nc")' ) TRIM( c_asmbkg ) 
    … …  
    97172            CALL iom_open( c_asmbkg, inum, ldwrt = .TRUE., kiolib = jprstlib) 
    98173            ! 
    99             IF( nitbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0 
    100                zdate = REAL( ndastp ) 
    101 #if defined key_zdftke 
    102                ! lk_zdftke=T :   Read turbulent kinetic energy ( en ) 
    103                IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...' 
    104                CALL tke_rst( nit000, 'READ' )               ! lk_zdftke=T :   Read turbulent kinetic energy ( en ) 
    105  
    106 #endif 
     174            ! 
     175            ! Write the information 
     176            IF ( ln_avgbkg ) THEN 
     177               IF( nitavgbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitavgbkg = 0 
     178                  zdate = REAL( ndastp ) 
     179#if defined key_zdftke 
     180                  ! lk_zdftke=T :   Read turbulent kinetic energy ( en ) 
     181                  IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...' 
     182                  CALL tke_rst( nit000, 'READ' )               ! lk_zdftke=T :   Read turbulent kinetic energy ( en ) 
     183 
     184#endif 
     185               ELSE 
     186                  zdate = REAL( ndastp ) 
     187               ENDIF 
     188               CALL iom_rstput( kt, nitavgbkg_r, inum, 'rdastp' , zdate   ) 
     189               CALL iom_rstput( kt, nitavgbkg_r, inum, 'un'     , un_tavg ) 
     190               CALL iom_rstput( kt, nitavgbkg_r, inum, 'vn'     , vn_tavg ) 
     191               CALL iom_rstput( kt, nitavgbkg_r, inum, 'tn'     , tn_tavg ) 
     192               CALL iom_rstput( kt, nitavgbkg_r, inum, 'sn'     , sn_tavg ) 
     193               CALL iom_rstput( kt, nitavgbkg_r, inum, 'sshn'   , sshn_tavg) 
     194#if defined key_zdftke 
     195               CALL iom_rstput( kt, nitavgbkg_r, inum, 'en'     , en_tavg ) 
     196#endif 
     197               CALL iom_rstput( kt, nitavgbkg_r, inum, 'avt'    , avt_tavg) 
     198               ! 
    107199            ELSE 
    108                zdate = REAL( ndastp ) 
     200               IF( nitbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0 
     201                  zdate = REAL( ndastp ) 
     202#if defined key_zdftke 
     203                  ! lk_zdftke=T :   Read turbulent kinetic energy ( en ) 
     204                  IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...' 
     205                  CALL tke_rst( nit000, 'READ' )               ! lk_zdftke=T :   Read turbulent kinetic energy ( en ) 
     206 
     207#endif 
     208               ELSE 
     209                  zdate = REAL( ndastp ) 
     210               ENDIF 
     211               CALL iom_rstput( kt, nitbkg_r, inum, 'rdastp' , zdate   ) 
     212               CALL iom_rstput( kt, nitbkg_r, inum, 'un'     , un                ) 
     213               CALL iom_rstput( kt, nitbkg_r, inum, 'vn'     , vn                ) 
     214               CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , tsn(:,:,:,jp_tem) ) 
     215               CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , tsn(:,:,:,jp_sal) ) 
     216               CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , sshn              ) 
     217#if defined key_zdftke 
     218               CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                ) 
     219#endif 
     220               CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               ) 
     221               ! 
    109222            ENDIF 
    110             ! 
    111             !                                      ! Write the information 
    112             CALL iom_rstput( kt, nitbkg_r, inum, 'rdastp' , zdate             ) 
    113             CALL iom_rstput( kt, nitbkg_r, inum, 'un'     , un                ) 
    114             CALL iom_rstput( kt, nitbkg_r, inum, 'vn'     , vn                ) 
    115             CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , tsn(:,:,:,jp_tem) ) 
    116             CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , tsn(:,:,:,jp_sal) ) 
    117             CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , sshn              ) 
    118 #if defined key_zdftke 
    119             CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                ) 
    120 #endif 
    121             CALL iom_rstput( kt, nitbkg_r, inum, 'gcx'    , gcx               ) 
    122             ! 
     223             
    123224            CALL iom_close( inum ) 
     225 
    124226         ENDIF 
    125227         ! 

  • branches/UKMO/CCRS_NUS_MC_STABLE/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r8058 r10177  
    4040#endif 
    4141   USE sbc_oce          ! Surface boundary condition variables. 
     42   USE zdfmxl, ONLY :  &   
     43   &  hmld_tref,       &    
     44#if defined key_karaml 
     45   &  hmld_kara,       & 
     46   &  ln_kara,         & 
     47#endif    
     48   &  hmld,            &  
     49   &  hmlp,            & 
     50   &  hmlpt 
     51#if defined key_bdy  
     52   USE bdy_oce, ONLY: bdytmask   
     53#endif   
    4254 
    4355   IMPLICIT NONE 
    … …  
    5769#endif 
    5870   LOGICAL, PUBLIC :: ln_bkgwri = .FALSE.      !: No output of the background state fields 
     71   LOGICAL, PUBLIC :: ln_avgbkg = .FALSE.      !: No output of the mean background state fields 
    5972   LOGICAL, PUBLIC :: ln_asmiau = .FALSE.      !: No applying forcing with an assimilation increment 
    6073   LOGICAL, PUBLIC :: ln_asmdin = .FALSE.      !: No direct initialization 
    … …  
    8093   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval  
    8194   INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval 
     95   INTEGER , PUBLIC ::   nitavgbkg   !: Number of timesteps to average assim bkg [0,nitavgbkg] 
    8296   !  
    8397   INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting 
    … …  
    87101   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment 
    88102   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc 
     103 
     104   INTEGER :: mld_choice        = 4   !: choice of mld criteria to use for physics assimilation 
     105                                      !: 1) hmld      - Turbocline/mixing depth                           [W points] 
     106                                      !: 2) hmlp      - Density criterion (0.01 kg/m^3 change from 10m)   [W points] 
     107                                      !: 3) hmld_kara - Kara MLD                                          [Interpolated] 
     108                                      !: 4) hmld_tref - Temperature criterion (0.2 K change from surface) [T points] 
     109 
    89110 
    90111   !! * Substitutions 
    … …  
    119140      INTEGER :: iitiaustr_date  ! Date YYYYMMDD of IAU interval start time step 
    120141      INTEGER :: iitiaufin_date  ! Date YYYYMMDD of IAU interval final time step 
     142      INTEGER :: isurfstat       ! Local integer for status of reading surft variable 
     143      INTEGER :: iitavgbkg_date  ! Date YYYYMMDD of end of assim bkg averaging period 
    121144      ! 
    122145      REAL(wp) :: znorm        ! Normalization factor for IAU weights 
    … …  
    127150      REAL(wp) :: zdate_inc    ! Time axis in increments file 
    128151      ! 
     152      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: &  
     153          &       t_bkginc_2d  ! file for reading in 2D   
     154      !                        ! temperature increments  
     155      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: &  
     156          &       z_mld     ! Mixed layer depth  
     157           
    129158      REAL(wp), POINTER, DIMENSION(:,:) ::   hdiv   ! 2D workspace 
    130       !! 
    131       NAMELIST/nam_asminc/ ln_bkgwri,                                      & 
     159      ! 
     160      LOGICAL :: lk_surft      ! Logical: T => Increments file contains surft variable  
     161                               !               so only apply surft increments. 
     162      !! 
     163      NAMELIST/nam_asminc/ ln_bkgwri, ln_avgbkg,                           & 
    132164         &                 ln_trainc, ln_dyninc, ln_sshinc,                & 
    133165         &                 ln_asmdin, ln_asmiau,                           & 
    134166         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   & 
    135          &                 ln_salfix, salfixmin, nn_divdmp 
     167         &                 ln_salfix, salfixmin, nn_divdmp, nitavgbkg, mld_choice 
    136168      !!---------------------------------------------------------------------- 
    137169 
    … …  
    139171      ! Read Namelist nam_asminc : assimilation increment interface 
    140172      !----------------------------------------------------------------------- 
     173 
     174      ! Set default values 
     175      ln_bkgwri = .FALSE. 
     176      ln_avgbkg = .FALSE. 
     177      ln_trainc = .FALSE. 
     178      ln_dyninc = .FALSE. 
     179      ln_sshinc = .FALSE. 
    141180      ln_seaiceinc = .FALSE. 
     181      ln_asmdin = .FALSE. 
     182      ln_asmiau = .TRUE. 
     183      ln_salfix = .FALSE. 
    142184      ln_temnofreeze = .FALSE. 
     185      salfixmin = -9999 
     186      nitbkg    = 0 
     187      nitdin    = 0       
     188      nitiaustr = 1 
     189      nitiaufin = 150 
     190      niaufn    = 0 
     191      nitavgbkg = 1 
    143192 
    144193      REWIND( numnam_ref )              ! Namelist nam_asminc in reference namelist : Assimilation increment 
    … …  
    158207         WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters' 
    159208         WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri 
     209         WRITE(numout,*) '      Logical switch for writing mean background state         ln_avgbkg = ', ln_avgbkg 
    160210         WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc 
    161211         WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc 
    … …  
    168218         WRITE(numout,*) '      Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr 
    169219         WRITE(numout,*) '      Timestep of end of IAU interval in [0,nitend-nit000-1]   nitiaufin = ', nitiaufin 
     220         WRITE(numout,*) '      Number of timesteps to average assim bkg [0,nitavgbkg]   nitavgbkg = ', nitavgbkg 
    170221         WRITE(numout,*) '      Type of IAU weighting function                           niaufn    = ', niaufn 
    171222         WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix 
    172223         WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin 
     224         WRITE(numout,*) '      Choice of MLD for physics assimilation                  mld_choice = ', mld_choice 
    173225      ENDIF 
    174226 
    … …  
    177229      nitiaustr_r = nitiaustr + nit000 - 1  ! Start of IAU interval referenced to nit000 
    178230      nitiaufin_r = nitiaufin + nit000 - 1  ! End of IAU interval referenced to nit000 
     231      nitavgbkg_r = nitavgbkg + nit000 - 1  ! Averaging period referenced to nit000 
    179232 
    180233      iiauper = nitiaufin_r - nitiaustr_r + 1  ! IAU interval length 
    … …  
    186239      CALL calc_date( nit000, nitiaustr_r, ndate0, iitiaustr_date )     ! IAU start time referenced to ndate0 
    187240      CALL calc_date( nit000, nitiaufin_r, ndate0, iitiaufin_date )     ! IAU end time referenced to ndate0 
     241      CALL calc_date( nit000, nitavgbkg_r, ndate0, iitavgbkg_date )     ! End of assim bkg averaging period referenced to ndate0 
    188242      ! 
    189243      IF(lwp) THEN 
    … …  
    197251         WRITE(numout,*) '       nitiaustr_r = ', nitiaustr_r 
    198252         WRITE(numout,*) '       nitiaufin_r = ', nitiaufin_r 
     253         WRITE(numout,*) '       nitavgbkg_r = ', nitavgbkg_r 
    199254         WRITE(numout,*) 
    200255         WRITE(numout,*) '   Dates referenced to current cycle:' 
    … …  
    206261         WRITE(numout,*) '       iitiaustr_date = ', iitiaustr_date 
    207262         WRITE(numout,*) '       iitiaufin_date = ', iitiaufin_date 
     263         WRITE(numout,*) '       iitavgbkg_date = ', iitavgbkg_date 
    208264      ENDIF 
    209265 
    … …  
    248304         & CALL ctl_stop( ' nitdin :',  & 
    249305         &                ' Background time step for Direct Initialization is outside', & 
     306         &                ' the cycle interval') 
     307 
     308      IF ( nitavgbkg_r > nitend ) & 
     309         & CALL ctl_stop( ' nitavgbkg_r :',  & 
     310         &                ' Assim bkg averaging period is outside', & 
    250311         &                ' the cycle interval') 
    251312 
    … …  
    327388      !-------------------------------------------------------------------- 
    328389 
    329       ALLOCATE( t_bkginc(jpi,jpj,jpk) ) 
    330       ALLOCATE( s_bkginc(jpi,jpj,jpk) ) 
    331       ALLOCATE( u_bkginc(jpi,jpj,jpk) ) 
    332       ALLOCATE( v_bkginc(jpi,jpj,jpk) ) 
    333       ALLOCATE( ssh_bkginc(jpi,jpj)   ) 
    334       ALLOCATE( seaice_bkginc(jpi,jpj)) 
     390      IF ( ln_trainc ) THEN 
     391         ALLOCATE( t_bkginc(jpi,jpj,jpk) ) 
     392         ALLOCATE( s_bkginc(jpi,jpj,jpk) ) 
     393         t_bkginc(:,:,:) = 0.0 
     394         s_bkginc(:,:,:) = 0.0 
     395      ENDIF 
     396      IF ( ln_dyninc ) THEN  
     397         ALLOCATE( u_bkginc(jpi,jpj,jpk) ) 
     398         ALLOCATE( v_bkginc(jpi,jpj,jpk) ) 
     399         u_bkginc(:,:,:) = 0.0 
     400         v_bkginc(:,:,:) = 0.0 
     401      ENDIF 
     402      IF ( ln_sshinc ) THEN 
     403         ALLOCATE( ssh_bkginc(jpi,jpj)   ) 
     404         ssh_bkginc(:,:) = 0.0 
     405      ENDIF 
     406      IF ( ln_seaiceinc ) THEN  
     407         ALLOCATE( seaice_bkginc(jpi,jpj)) 
     408         seaice_bkginc(:,:) = 0.0 
     409      ENDIF 
    335410#if defined key_asminc 
    336411      ALLOCATE( ssh_iau(jpi,jpj)      ) 
    337 #endif 
    338       t_bkginc(:,:,:) = 0.0 
    339       s_bkginc(:,:,:) = 0.0 
    340       u_bkginc(:,:,:) = 0.0 
    341       v_bkginc(:,:,:) = 0.0 
    342       ssh_bkginc(:,:) = 0.0 
    343       seaice_bkginc(:,:) = 0.0 
    344 #if defined key_asminc 
    345412      ssh_iau(:,:)    = 0.0 
    346413#endif 
    … …  
    378445 
    379446         IF ( ln_trainc ) THEN    
    380             CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 ) 
    381             CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 ) 
    382             ! Apply the masks 
    383             t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:) 
    384             s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:) 
    385             ! Set missing increments to 0.0 rather than 1e+20 
    386             ! to allow for differences in masks 
    387             WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0 
    388             WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0 
     447 
     448            !Test if the increments file contains the surft variable. 
     449            isurfstat = iom_varid( inum, 'bckinsurft', ldstop = .FALSE. ) 
     450            IF ( isurfstat == -1 ) THEN 
     451               lk_surft = .FALSE. 
     452            ELSE 
     453               lk_surft = .TRUE. 
     454               CALL ctl_warn( ' Applying 2D temperature increment to bottom of ML: ', & 
     455            &                 ' bckinsurft found in increments file.' ) 
     456            ENDIF              
     457 
     458            IF (lk_surft) THEN  
     459                 
     460               ALLOCATE(z_mld(jpi,jpj))  
     461               SELECT CASE(mld_choice)  
     462               CASE(1)  
     463                  z_mld = hmld  
     464               CASE(2)  
     465                  z_mld = hmlp  
     466               CASE(3)  
     467#if defined key_karaml 
     468                  IF ( ln_kara ) THEN 
     469                     z_mld = hmld_kara 
     470                  ELSE 
     471                     CALL ctl_stop("Kara mixed layer not calculated as ln_kara=.false.") 
     472                  ENDIF 
     473#else 
     474                  CALL ctl_stop("Kara mixed layer not defined in current version of NEMO")  ! JW: Safety feature, should be removed 
     475                                                                                            ! once the Kara mixed layer is available  
     476#endif 
     477               CASE(4)  
     478                  z_mld = hmld_tref  
     479               END SELECT        
     480                       
     481               ALLOCATE( t_bkginc_2d(jpi,jpj) )  
     482               CALL iom_get( inum, jpdom_autoglo, 'bckinsurft', t_bkginc_2d, 1)  
     483#if defined key_bdy                 
     484               DO jk = 1,jpkm1  
     485                  WHERE( z_mld(:,:) > fsdepw(:,:,jk) )  
     486                     t_bkginc(:,:,jk) = t_bkginc_2d(:,:) * 0.5 * & 
     487                     &       ( 1 + cos( (fsdept(:,:,jk)/z_mld(:,:) ) * rpi ) ) 
     488                      
     489                     t_bkginc(:,:,jk) = t_bkginc(:,:,jk) * bdytmask(:,:)  
     490                  ELSEWHERE  
     491                     t_bkginc(:,:,jk) = 0.  
     492                  ENDWHERE  
     493               ENDDO  
     494#else  
     495               t_bkginc(:,:,:) = 0.  
     496#endif                 
     497               s_bkginc(:,:,:) = 0.  
     498                 
     499               DEALLOCATE(z_mld, t_bkginc_2d)  
     500             
     501            ELSE  
     502                
     503               CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 ) 
     504               CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 ) 
     505               ! Apply the masks 
     506               t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:) 
     507               s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:) 
     508               ! Set missing increments to 0.0 rather than 1e+20 
     509               ! to allow for differences in masks 
     510               WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0 
     511               WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0 
     512          
     513            ENDIF 
     514 
    389515         ENDIF 
    390516 
    … …  
    8921018            ENDIF 
    8931019 
     1020#if defined key_asminc 
     1021         ELSE 
     1022            ssh_iau(:,:) = 0._wp 
     1023#endif 
    8941024         ENDIF 
    8951025 

  • branches/UKMO/CCRS_NUS_MC_STABLE/NEMOGCM/NEMO/OPA_SRC/ASM/asmpar.F90

    r8058 r10177  
    2929   INTEGER, PUBLIC :: nitiaufin_r   !: IAU final time step referenced to nit000 
    3030   INTEGER, PUBLIC :: nittrjfrq     !: Frequency of trajectory output for 4D-VAR 
     31   INTEGER, PUBLIC :: nitavgbkg_r   !: Averaging period for assim bkg referenced to nit000 
    3132 
    3233   !!---------------------------------------------------------------------- 

  • branches/UKMO/CCRS_NUS_MC_STABLE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r8059 r10177  
    4646   USE agrif_opa_interp ! agrif 
    4747#endif 
    48 #if defined key_asminc    
    49    USE asminc          ! Assimilation increment 
    50 #endif 
    5148 
    5249   IMPLICIT NONE 
    … …  
    462459                &                        + fwfisf(:,:) + fwfisf_b(:,:)                     ) 
    463460      ENDIF 
    464 #if defined key_asminc 
    465       !                                         ! Include the IAU weighted SSH increment 
    466       IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 
    467          zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:) 
    468       ENDIF 
    469 #endif 
    470461      !                                   !* Fill boundary data arrays for AGRIF 
    471462      !                                   ! ------------------------------------ 

  • branches/UKMO/CCRS_NUS_MC_STABLE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r8058 r10177  
    7474      INTEGER, INTENT(in) ::   kt                      ! time step 
    7575      !  
    76       INTEGER             ::   jk                      ! dummy loop indice 
     76      INTEGER             ::   ji, jj, jk              ! dummy loop indices 
    7777      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars 
    7878      !!---------------------------------------------------------------------- 
    … …  
    9494      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog) 
    9595      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt 
     96 
     97#if defined key_asminc 
     98      !                                                ! Include the IAU weighted SSH increment 
     99      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 
     100         CALL ssh_asm_inc( kt ) 
     101#if defined key_vvl 
     102! Don't directly adjust ssh but change hdivn at all levels instead 
     103! In trasbc also add in the heat and salt content associated with these changes at each level   
     104        DO jk = 1, jpkm1                                  
     105                 hdivn(:,:,jk) = hdivn(:,:,jk) - ( ssh_iau(:,:) / ( ht_0(:,:) + 1.0 - ssmask(:,:) ) ) * ( e3t_0(:,:,jk) / fse3t_n(:,:,jk) ) * tmask(:,:,jk)  
     106        END DO 
     107        CALL lbc_lnk( hdivn, 'T', 1. ) ! Not sure that's necessary 
     108      ENDIF 
     109#endif 
     110#endif 
    96111 
    97112      !                                           !------------------------------! 
    … …  
    123138#endif 
    124139 
    125 #if defined key_asminc 
    126       !                                                ! Include the IAU weighted SSH increment 
    127       IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 
    128          CALL ssh_asm_inc( kt ) 
    129          ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:) 
    130       ENDIF 
    131 #endif 
    132  
    133140      !                                           !------------------------------! 
    134141      !                                           !           outputs            ! 

  • branches/UKMO/CCRS_NUS_MC_STABLE/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

    r8059 r10177  
    121121      IF( lk_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries 
    122122#endif 
     123  
     124      ! set time step size (Euler/Leapfrog) 
     125      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dtra(:) =     rdttra(:)      ! at nit000             (Euler) 
     126      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dtra(:) = 2._wp* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog) 
     127      ENDIF 
    123128 
    124129#if ( ! defined key_lim3 && ! defined key_lim2 && ! key_cice ) 
    … …  
    151156#endif 
    152157 
    153       ! set time step size (Euler/Leapfrog) 
    154       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dtra(:) =     rdttra(:)      ! at nit000             (Euler) 
    155       ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dtra(:) = 2._wp* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog) 
    156       ENDIF 
    157  
    158158      ! trends computation initialisation 
    159159      IF( l_trdtra )   THEN                    ! store now fields before applying the Asselin filter 

  • branches/UKMO/CCRS_NUS_MC_STABLE/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r8059 r10177  
    3434   USE timing          ! Timing 
    3535   USE eosbn2 
     36#if defined key_asminc    
     37   USE asminc          ! Assimilation increment 
     38#endif 
    3639 
    3740   IMPLICIT NONE 
    … …  
    279282         END DO   
    280283      ENDIF 
     284 
     285#if defined key_asminc 
     286! WARNING: THIS MAY WELL NOT BE REQUIRED - WE DON'T WANT TO CHANGE T&S BUT THIS MAY COMPENSATE ANOTHER TERM... 
     287! Rate of change in e3t for each level is ssh_iau*e3t_0/ht_0 
     288! Contribution to tsa should be rate of change in level / per m of ocean? (hence the division by fse3t_n) 
     289      IF( ln_sshinc ) THEN         ! input of heat and salt due to assimilation 
     290         DO jj = 2, jpj  
     291            DO ji = fs_2, fs_jpim1 
     292               zdep = ssh_iau(ji,jj) / ( ht_0(ji,jj) + 1.0 - ssmask(ji, jj) ) 
     293               DO jk = 1, jpkm1 
     294                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   & 
     295                                        &            + tsn(ji,jj,jk,jp_tem) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) ) 
     296                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   & 
     297                                        &            + tsn(ji,jj,jk,jp_sal) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) ) 
     298               END DO 
     299            END DO   
     300         END DO   
     301      ENDIF 
     302#endif 
    281303  
    282304      IF( l_trdtra )   THEN                      ! send trends for further diagnostics 

  • branches/UKMO/CCRS_NUS_MC_STABLE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

    r8393 r10177  
    2727   PRIVATE 
    2828 
     29   PUBLIC   zdf_mxl_tref  ! called by asminc.F90 
    2930   PUBLIC   zdf_mxl       ! called by step.F90 
    3031 
     32   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld_tref  !: mixed layer depth at t-points - temperature criterion [m] 
    3133   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln    !: number of level in the mixed layer (used by TOP) 
    3234   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld    !: mixing layer depth (turbocline)      [m] 
    … …  
    7880        &          ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc ) 
    7981         ! 
     82         ALLOCATE(hmld_tref(jpi,jpj)) 
    8083         IF( lk_mpp             )   CALL mpp_sum ( zdf_mxl_alloc ) 
    8184         IF( zdf_mxl_alloc /= 0 )   CALL ctl_warn('zdf_mxl_alloc: failed to allocate arrays.') 
    … …  
    8487   END FUNCTION zdf_mxl_alloc 
    8588 
     89   SUBROUTINE zdf_mxl_tref()  
     90      !!----------------------------------------------------------------------  
     91      !!                  ***  ROUTINE zdf_mxl_tref  ***  
     92      !!                     
     93      !! ** Purpose :   Compute the mixed layer depth with temperature criteria.  
     94      !!  
     95      !! ** Method  :   The temperature-defined mixed layer depth is required  
     96      !!                   when assimilating SST in a 2D analysis.   
     97      !!  
     98      !! ** Action  :   hmld_tref  
     99      !!----------------------------------------------------------------------  
     100      !  
     101      INTEGER  ::   ji, jj, jk   ! dummy loop indices  
     102      REAL(wp) ::   t_ref               ! Reference temperature    
     103      REAL(wp) ::   temp_c = 0.2        ! temperature criterion for mixed layer depth    
     104      !!----------------------------------------------------------------------  
     105      !  
     106      ! Initialise array  
     107      IF( zdf_mxl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl_tref : unable to allocate arrays' )  
     108        
     109      !For the AMM model assimiation uses a temperature based mixed layer depth    
     110      !This is defined here    
     111      DO jj = 1, jpj    
     112         DO ji = 1, jpi    
     113           hmld_tref(ji,jj)=fsdept(ji,jj,1  )     
     114           IF(ssmask(ji,jj) > 0.)THEN    
     115             t_ref=tsn(ji,jj,1,jp_tem)   
     116             DO jk=2,jpk    
     117               IF(ssmask(ji,jj)==0.)THEN    
     118                  hmld_tref(ji,jj)=fsdept(ji,jj,jk )    
     119                  EXIT    
     120               ELSEIF( ABS(tsn(ji,jj,jk,jp_tem)-t_ref) < temp_c)THEN    
     121                  hmld_tref(ji,jj)=fsdept(ji,jj,jk )    
     122               ELSE    
     123                  EXIT    
     124               ENDIF    
     125             ENDDO    
     126           ENDIF    
     127         ENDDO    
     128      ENDDO  
     129    
     130   END SUBROUTINE zdf_mxl_tref  
    86131 
    87132   SUBROUTINE zdf_mxl( kt ) 

  • branches/UKMO/CCRS_NUS_MC_STABLE/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r8561 r10177  
    480480 
    481481      !                                     ! Assimilation increments 
    482       IF( lk_asminc     )   CALL asm_inc_init   ! Initialize assimilation increments 
     482      IF( lk_asminc ) THEN  
     483#if defined key_shelf  
     484         CALL  zdf_mxl_tref()     ! Initialization of hmld_tref 
     485#endif  
     486         CALL asm_inc_init     ! Initialize assimilation increments  
     487      ENDIF 
     488 
    483489      IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler 
    484490                            CALL dia_tmb_init  ! TMB outputs 

  • branches/UKMO/CCRS_NUS_MC_STABLE/NEMOGCM/NEMO/OPA_SRC/step.F90

    r8561 r10177  
    336336      ! Dynamics                                    (tsa used as workspace) 
    337337      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     338 
     339      IF( ln_bkgwri )        CALL asm_bkg_wri( kstp )     ! output background fields 
     340 
    338341      IF( lk_dynspg_ts   )  THEN 
    339342                                                             ! revert to previously computed momentum tendencies 
    … …  
    354357        IF(  lk_asminc .AND. ln_asmiau .AND. & 
    355358           & ln_dyninc      )  CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment 
    356         IF( ln_bkgwri )        CALL asm_bkg_wri( kstp )     ! output background fields 
    357359        IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified) 
    358360        IF( lk_bdy          )  CALL bdy_dyn3d_dmp(kstp )    ! bdy damping trends 
    … …  
    373375                               CALL ssh_swp( kstp )         ! swap of sea surface height 
    374376      IF( lk_vvl           )   CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors 
    375       ! 
    376377 
    377378#if defined key_agrif 
    … …  
    403404         IF( lwm.AND.numoni /= -1 ) CALL FLUSH    ( numoni )     ! flush output namelist ice 
    404405      ENDIF 
    405       IF( lrst_oce         )   CALL rst_write( kstp )       ! write output ocean restart file 
     406      IF( lrst_oce         )   CALL rst_write    ( kstp )   ! write output ocean restart file 
    406407 
    407408      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
Note: See TracChangeset for help on using the changeset viewer.