Ignore:
Timestamp:
10/05/16 14:25:40 (8 years ago)
Author:
vancop
Message:

Ludivine source files

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_bio_sms.f

    r4 r20  
    8282!-----------------------------------------------------------------------------! 
    8383! 
     84      consN = 0 
     85 
     86      DO layer = layer_00, nlay_bio  
     87          
     88         trucN = cbu_i_bio(jn_din,layer) + cbu_i_bio(jn_eon,layer) 
     89     &              + cbu_i_bio(jn_aon,layer) + consN 
     90         consN    = trucN         
     91 
     92      END DO 
     93 
     94      WRITE(numout,*) ' Conservation BIO N :', consN 
     95 
    8496 
    8597      ! Reference of the upper layer 
     
    246258      END DO 
    247259 
    248       ! Remineralization  
    249       IF ( ln_rem ) THEN 
    250  
    251          SELECT CASE ( nn_bio_opt )  
    252  
    253          CASE(0) ! NP 
    254  
    255             rem_bio(layer_00:nlay_bio) = frem_bio *  
    256      &                                   lys_bio(layer_00:nlay_bio) 
    257  
    258          CASE(1) ! NPDr 
    259  
    260             IF ( nn_rem .GE. 1 )   ! no T-dependence 
    261      &         rem_bio(layer_00:nlay_bio) = krem_bio *  
    262      &                          cbu_i_bio(jn_eoc,layer_00:nlay_bio) 
    263  
    264           
    265             IF ( nn_rem .EQ. 2 )   ! T-dependence 
    266      &         rem_bio(layer_00:nlay_bio) = rem_bio(layer_00:nlay_bio) * 
    267      &                       EXP( rg_bac * tc_bio(layer_00:nlay_bio) ) 
    268  
    269          END SELECT 
    270  
    271       ENDIF 
    272  
     260      DO layer = layer_00, nlay_bio 
     261 
     262         IF ( ln_rem ) THEN 
     263 
     264            SELECT CASE ( nn_bio_opt )  
     265 
     266               CASE(0) ! NP 
     267 
     268                  rem_bio(layer) = frem_bio *  
     269     &                             lys_bio(layer) 
     270 
     271               CASE(1) ! NPD 
     272  
     273                  IF ( nn_rem .EQ. 1 )   ! no T-dependence 
     274     &               rem_bio(layer) = krem_bio *  
     275     &                                cbu_i_bio(jn_eoc,layer) 
     276 
     277                  IF ( nn_rem .EQ. 2 )   ! T-dependence 
     278     &               rem_bio(layer) = krem_bio * cbu_i_bio(jn_eoc,layer) 
     279     &                                * EXP( rg_bac * tc_bio(layer) ) 
     280      
     281            END SELECT 
     282 
     283         ENDIF 
     284 
     285      END DO 
    273286 
    274287      !------------------------------------------------------------------------ 
    275288      ! 2.3 Diatom Synthesis 
    276289      !------------------------------------------------------------------------ 
     290      zeps = 0.000000001 
    277291 
    278292      DO layer = layer_00, nlay_bio 
    279     
     293 
    280294         IF ( ln_syn ) THEN 
    281295            ! raw synthesis (in s-1) 
     
    287301            zsyn1 = z_syn * cbu_i_bio(jn_aoc,layer) ! normal carbon production rate (mmolC/m3/s) 
    288302            zsyn2 = zsyn1;    zsyn3 = zsyn1;    zsyn4 = zsyn1 
     303 
    289304    
    290305            IF ( ln_lim_dsi )                  ! prevent exhaustion of Si 
    291306     &         zsyn2 = cbu_i_bio(jn_dsi,layer) / ddtb / si_c - zeps  
    292307   
    293             IF ( ln_lim_no3 )                  ! prevent exhaustion of N 
    294      &         zsyn3 = cbu_i_bio(jn_din,layer) / ddtb / no3_c - zeps 
    295    
    296             IF ( ln_lim_po4 )                  ! prevent exhaustion of P 
    297      &         zsyn4 = cbu_i_bio(jn_dip,layer) / ddtb / po4_c - zeps 
    298     
    299             ! nutrient stock-limited PP rate 
     308            IF (ln_decoupNC) THEN ! C aussi limité en stock par DIN 
     309               zC  = MAX (cbu_i_bio(jn_aoc,layer) , zeps ) 
     310               zN  = MAX (cbu_i_bio(jn_aon,layer) , zeps ) 
     311               N_C_alg(layer)  = zN / zC 
     312!              N_C_alg(layer) = 0.14  
     313               zsyn3 = cbu_i_bio(jn_din,layer) / (ddtb*N_C_alg(layer)) 
     314     &                 - zeps 
     315            ELSE  
     316               zsyn3 = cbu_i_bio(jn_din,layer) / (ddtb * no3_c) - zeps 
     317            ENDIF 
     318 
     319 
     320            IF (ln_decoupPC) THEN ! C aussi limité en stock par DIP 
     321               zC  = MAX (cbu_i_bio(jn_aoc,layer) , zeps ) 
     322               zP  = MAX (cbu_i_bio(jn_aop,layer) , zeps ) 
     323               P_C_alg(layer)  = zP / zC 
     324!              P_C_alg(layer) = 0.0094  
     325               zsyn4 = cbu_i_bio(jn_dip,layer)/(ddtb * P_C_alg(layer)) 
     326     &                 - zeps 
     327            ELSE 
     328               zsyn4 = cbu_i_bio(jn_dip,layer) / ddtb / po4_c - zeps 
     329            ENDIF 
     330 
     331 
     332             ! nutrient stock-limited PP rate 
    300333            syn_bio(layer) = MIN( zsyn1, zsyn2, zsyn3, zsyn4 )  
    301334 
     
    312345    
    313346 
    314       END DO ! layer 
    315  
    316       !------------------------------------------------------------------------ 
    317       ! 2.4 Update reservoirs 
    318       !------------------------------------------------------------------------ 
    319  
    320       DO layer = layer_00, nlay_bio 
     347         ! ------------------------- 
     348         ! -- Limitations en stocks 
     349         ! ------------------------ 
     350 
     351         IF (ln_decoupNC) THEN 
     352             zzdin = syn_bio(layer) / (((1/N_C_alg(layer)) * 
     353     &               (cbu_i_bio(jn_din,layer)/ddtb)) - zeps) 
     354 
     355         ELSE  
     356             zzdin = syn_bio(layer) / (((1/no3_c) * 
     357     &               (cbu_i_bio(jn_din,layer)/ddtb)) - zeps) 
     358         ENDIF 
     359 
     360 
     361         IF (ln_decoupPC) THEN 
     362             zzpo4 = syn_bio(layer) / (((1/P_C_alg(layer)) * 
     363     &               (cbu_i_bio(jn_dip,layer)/ddtb)) - zeps) 
     364 
     365         ELSE 
     366             zzpo4 = syn_bio(layer) / (((1/po4_c) * 
     367     &               (cbu_i_bio(jn_dip,layer)/ddtb)) - zeps) 
     368         ENDIF 
     369 
     370 
     371 
     372         zzsi  = syn_bio(layer) / (((1/si_c) * 
     373     &           (cbu_i_bio(jn_dsi,layer)/ddtb)) - zeps) 
     374 
     375 
     376         IF ( zzdin .EQ. 1 ) THEN 
     377            lim_din_stock(layer) = 1   !DIN est limitant en stocks 
     378         ELSE 
     379            lim_din_stock(layer) = 0   !DIN PAS limitant en stocks 
     380         ENDIF 
     381 
     382         IF ( zzpo4 .EQ. 1 ) THEN 
     383            lim_dip_stock(layer) = 1 
     384         ELSE 
     385            lim_dip_stock(layer) = 0 
     386         ENDIF 
     387 
     388         IF ( zzsi .EQ. 1) THEN 
     389            lim_dsi_stock(layer) = 1 
     390         ELSE 
     391            lim_dsi_stock(layer) = 0 
     392         ENDIF 
     393 
     394 
     395      !------------------------------------------------------------------------ 
     396      ! N cycle  
     397      !------------------------------------------------------------------------ 
     398 
     399 
     400         IF ( ln_decoupNC ) THEN 
     401 
     402            !--- Nitrogen carbon algal ratio 
     403            zC    = MAX( cbu_i_bio(jn_aoc,layer) , zeps ) 
     404            zN    = MAX( cbu_i_bio(jn_aon,layer) , zeps ) 
     405            N_C_alg(layer)  = zN / zC 
     406!           zN_C  = MIN( cbu_i_bio(jn_aon,layer) / zC , N_C_max ) 
     407!           zN_C  = MAX( N_C_min , zN_C ) 
     408!           N_C_alg(layer) = zN_C  ! => marche pas 
     409!            N_C_alg(layer) = 0.14 
     410 
     411            !--- Nitrogen carbon detritic ratio 
     412            zC    = MAX( cbu_i_bio(jn_eoc,layer), zeps ) 
     413            zNd   = MAX( cbu_i_bio(jn_eon,layer), zeps )  
     414            N_C_det(layer) = zNd / zC 
     415!            N_C_det(layer) = 0.14 
     416 
     417            !--- Process rates 
     418            syn_N(layer) = N_C_alg(layer) * ksyn_N * syn_bio(layer) 
     419            rsp_N(layer) = N_C_alg(layer) * krsp_N * rsp_bio(layer) 
     420            lys_N(layer) = N_C_alg(layer) * klys_N * lys_bio(layer) 
     421            rem_N(layer) = N_C_det(layer) * krem_N * rem_bio(layer) 
     422 
     423 
     424         ENDIF ! ln_decoupNC 
     425 
     426      !------------------------------------------------------------------------ 
     427      ! P cycle  
     428      !------------------------------------------------------------------------ 
     429 
     430 
     431         IF ( ln_decoupPC ) THEN 
     432 
     433            !--- P carbon algal ratio 
     434            zC    = MAX( cbu_i_bio(jn_aoc,layer) , zeps ) 
     435            zP    = MAX( cbu_i_bio(jn_aop,layer) , zeps ) 
     436            P_C_alg(layer)  = zP / zC 
     437 
     438            !--- P carbon detritic ratio 
     439            zC    = MAX( cbu_i_bio(jn_eoc,layer), zeps ) 
     440            zPd   = MAX( cbu_i_bio(jn_eop,layer), zeps ) 
     441            P_C_det(layer) = zPd / zC 
     442 
     443            !--- Process rates 
     444            syn_P(layer) = P_C_alg(layer) * ksyn_P * syn_bio(layer) 
     445            rsp_P(layer) = P_C_alg(layer) * krsp_P * rsp_bio(layer) 
     446            lys_P(layer) = P_C_alg(layer) * klys_P * lys_bio(layer) 
     447            rem_P(layer) = P_C_det(layer) * krem_P * rem_bio(layer) 
     448 
     449         ENDIF ! ln_decoupPC 
     450 
     451 
     452         zPa   = MAX( cbu_i_bio(jn_aop,layer) , zeps ) 
     453         zNt   = MAX( cbu_i_bio(jn_aon,layer) , zeps ) 
     454         N_P_alg(layer)  = zNt / zPa 
     455         zPe   = MAX( cbu_i_bio(jn_eop,layer) , zeps ) 
     456         zNtd  = MAX( cbu_i_bio(jn_eon,layer) , zeps ) 
     457         N_P_det(layer)  = zNtd / zPe 
     458 
     459 
     460 
     461      !------------------------------------------------------------------------ 
     462      ! 2.4 Update reservoirs (including N cycle reservoirs) 
     463      !------------------------------------------------------------------------ 
     464 
     465        ! ------------------------------------------------------- 
     466        ! -- 1) Increments 
     467        ! ------------------------------------------------------ 
     468 
    321469 
    322470         ! DSi uptake 
     
    325473         zd_dsi2 = - cbu_i_bio(jn_dsi,layer)        ! Maximum DSi uptake 
    326474         zd_dsi  = MAX ( zd_dsi1, zd_dsi2 ) 
    327  
    328          ! NO3 uptake 
    329          zd_no31 = no3_c * ( rem_bio(layer) + 
    330      &             rsp_bio(layer) - syn_bio(layer) ) * ddtb 
    331          zd_no32 = - cbu_i_bio(jn_din,layer)        ! Max NO3 uptake 
    332          zd_no3  = MAX ( zd_no31, zd_no32 ) 
    333  
    334          ! update PO4 
    335          zd_po41 = po4_c * ( rem_bio(layer) + 
    336      &             rsp_bio(layer) - syn_bio(layer) ) * ddtb 
    337          zd_po42 = - cbu_i_bio(jn_dip,layer)        ! Max PO4 uptake 
    338          zd_po4  = MAX ( zd_po41, zd_po42 ) 
    339475 
    340476         ! update algae 
     
    350486            zd_eoc2 = - cbu_i_bio(jn_eoc,layer)     ! Max detrital carbon loss 
    351487            zd_eoc  = MAX( zd_eoc1, zd_eoc2 ) 
    352          ENDIF 
     488         ENDIF ! nn_bio_opt 
     489 
     490 
     491         IF (ln_decoupNC) THEN 
     492 
     493            !--- Increments 
     494            zd_din1 = ( rsp_N(layer) - syn_N(layer) + 
     495     &                  rem_N(layer) ) * ddtb 
     496            zd_din2 = - cbu_i_bio(jn_din,layer) !PEU PAS EPUISER PLUS 
     497                                                !QUE CE QUIL Y A !!!! 
     498            zd_din  = MAX ( zd_din1, zd_din2 )  ! interessant dans le cas 
     499                                                ! ou zd_din1 est négatif 
     500                                                ! zd_din1 peut pas être 
     501                                                ! plus petit que zd_din2  
     502 
     503            zd_aon1 = ( syn_N(layer) - rsp_N(layer) - lys_N(layer) ) * 
     504     &                ddtb 
     505            zd_aon2 = - cbu_i_bio(jn_aon,layer) 
     506            zd_aon  = MAX ( zd_aon1, zd_aon2 ) 
     507 
     508            zd_eon1 = ( lys_N(layer) - rem_N(layer) ) * ddtb 
     509            zd_eon2 = - cbu_i_bio(jn_eon,layer) 
     510            zd_eon  = MAX ( zd_eon1, zd_eon2 ) 
     511 
     512            !--- Correct rates to prevent leaks 
     513            IF ( zd_din .EQ. zd_din2 ) 
     514     &         syn_N(layer) = MIN( syn_N(layer), -zd_din2/ddtb) 
     515 
     516            IF ( zd_aon .EQ. zd_aon2 ) 
     517     &         lys_N(layer) = MIN( lys_N(layer), -zd_aon2/ddtb) 
     518 
     519            IF ( zd_eon .EQ. zd_eon2 ) 
     520     &         rem_N(layer) = MIN( rem_N(layer), -zd_eon2/ddtb) 
     521 
     522 
     523         ELSE 
     524            zd_din1 = no3_c * ( rem_bio(layer) + 
     525     &                rsp_bio(layer) - syn_bio(layer) ) * ddtb 
     526            zd_din2 = - cbu_i_bio(jn_din,layer)        ! Max NO3 uptake 
     527            zd_din  = MAX ( zd_din1, zd_din2 )      
     528              
     529         ENDIF ! ln_decoupNC 
     530 
     531!!!!!!!!!!!!!!!!!! 
     532 
     533         IF (ln_decoupPC) THEN 
     534 
     535            !--- Increments 
     536            zd_dip1 = ( rsp_P(layer) - syn_P(layer) + 
     537     &                  rem_P(layer) ) * ddtb 
     538            zd_dip2 = - cbu_i_bio(jn_dip,layer) !PEU PAS EPUISER PLUS 
     539                                                !QUE CE QUIL Y A !!!! 
     540            zd_dip  = MAX ( zd_dip1, zd_dip2 )  ! interessant dans le 
     541                                                ! ou zd_din1 est négatif 
     542                                                ! zd_din1 peut pas être 
     543                                                ! plus petit que zd_din2  
     544 
     545            zd_aop1 = ( syn_P(layer) - rsp_P(layer) - lys_P(layer) ) * 
     546     &                ddtb 
     547            zd_aop2 = - cbu_i_bio(jn_aop,layer) 
     548            zd_aop  = MAX ( zd_aop1, zd_aop2 ) 
     549 
     550            zd_eop1 = ( lys_P(layer) - rem_P(layer) ) * ddtb 
     551            zd_eop2 = - cbu_i_bio(jn_eop,layer) 
     552            zd_eop  = MAX ( zd_eop1, zd_eop2 ) 
     553 
     554            !--- Correct rates to prevent leaks 
     555            IF ( zd_dip .EQ. zd_dip2 ) 
     556     &         syn_P(layer) = MIN( syn_P(layer), -zd_dip2/ddtb) 
     557 
     558            IF ( zd_aop .EQ. zd_aop2 ) 
     559     &         lys_P(layer) = MIN( lys_P(layer), -zd_aop2/ddtb) 
     560 
     561            IF ( zd_eop .EQ. zd_eop2 ) 
     562     &         rem_P(layer) = MIN( rem_P(layer), -zd_eop2/ddtb) 
     563 
     564 
     565         ELSE 
     566            zd_dip1 = po4_c * ( rem_bio(layer) + 
     567     &                rsp_bio(layer) - syn_bio(layer) ) * ddtb 
     568            zd_dip2 = - cbu_i_bio(jn_dip,layer)        ! Max DIP uptake 
     569            zd_dip  = MAX ( zd_dip1, zd_dip2 ) 
     570 
     571         ENDIF ! ln_decoupPC 
     572 
     573 
     574        ! ------------------------------------------------------- 
     575        ! -- 2) Add Increments 
     576        ! ------------------------------------------------------ 
     577 
     578         IF (ln_decoupNC) THEN 
     579            cbu_i_bio(jn_din,layer) = cbu_i_bio(jn_din,layer) + zd_din 
     580            cbu_i_bio(jn_aon,layer) = cbu_i_bio(jn_aon,layer) + zd_aon 
     581            cbu_i_bio(jn_eon,layer) = cbu_i_bio(jn_eon,layer) + zd_eon 
     582         ELSE 
     583            cbu_i_bio(jn_din,layer) = cbu_i_bio(jn_din,layer) + zd_din 
     584         ENDIF 
     585 
     586         IF (ln_decoupPC) THEN 
     587            cbu_i_bio(jn_dip,layer) = cbu_i_bio(jn_dip,layer) + zd_dip 
     588            cbu_i_bio(jn_aop,layer) = cbu_i_bio(jn_aop,layer) + zd_aop 
     589            cbu_i_bio(jn_eop,layer) = cbu_i_bio(jn_eop,layer) + zd_eop 
     590         ELSE 
     591            cbu_i_bio(jn_dip,layer) = cbu_i_bio(jn_dip,layer) + zd_dip 
     592         ENDIF 
     593 
    353594 
    354595         cbu_i_bio(jn_dsi,layer) = cbu_i_bio(jn_dsi,layer) + zd_dsi 
    355          cbu_i_bio(jn_din,layer) = cbu_i_bio(jn_din,layer) + zd_no3 
    356          cbu_i_bio(jn_dip,layer) = cbu_i_bio(jn_dip,layer) + zd_po4 
    357596         cbu_i_bio(jn_aoc,layer) = cbu_i_bio(jn_aoc,layer) + zd_daf 
    358597 
    359598         IF ( nn_bio_opt .GE. 1 )  
    360599     &      cbu_i_bio(jn_eoc,layer) = cbu_i_bio(jn_eoc,layer) + zd_eoc 
     600 
    361601 
    362602         IF ( ln_write_bio ) THEN 
     
    375615            WRITE(numout,*) ' lim_tem   : ', lim_tem(layer) 
    376616            WRITE(numout,*) ' lim_sal   : ', lim_sal(layer) 
     617            WRITE(numout,*) ' lim_din_stock : ', lim_din_stock(layer) 
     618            WRITE(numout,*) ' lim_dip_stock : ', lim_dip_stock(layer) 
     619            WRITE(numout,*) ' lim_dsi_stock : ', lim_dsi_stock(layer) 
    377620            WRITE(numout,*) 
    378621            WRITE(numout,*) ' tc_bio    : ', tc_bio(layer)   
     
    381624            WRITE(numout,*) 
    382625            WRITE(numout,*) ' zd_dsi    : ', zd_dsi 
    383             WRITE(numout,*) ' zd_no3    : ', zd_no3 
    384626            WRITE(numout,*) ' zd_po4    : ', zd_po4 
    385627            WRITE(numout,*) ' zd_daf    : ', zd_daf 
    386628            WRITE(numout,*) ' zd_eoc    : ', zd_eoc 
     629            WRITE(numout,*) ' zd_din    : ', zd_din 
    387630            WRITE(numout,*) 
    388631            WRITE(numout,*) ' dSi       : ', cbu_i_bio(jn_dsi,layer) 
    389             WRITE(numout,*) ' dIN       : ', cbu_i_bio(jn_din,layer) 
     632            WRITE(numout,*) ' dIN       : ', cbu_i_bio(jn_din,layer)  
    390633            WRITE(numout,*) ' dIP       : ', cbu_i_bio(jn_dip,layer) 
    391634            WRITE(numout,*) ' AoC       : ', cbu_i_bio(jn_aoc,layer) 
    392635            WRITE(numout,*) ' eoC       : ', cbu_i_bio(jn_eoc,layer) 
    393             WRITE(numout,*) 
     636            WRITE(numout,*) ' krem_N    : ', krem_N 
     637 
     638         
     639            IF (ln_decoupNC) THEN  
     640               WRITE(numout,*) ' syn_N    : ', syn_N(layer) 
     641               WRITE(numout,*) ' rsp_N    : ', rsp_N(layer) 
     642               WRITE(numout,*) ' lys_N    : ', lys_N(layer) 
     643               WRITE(numout,*) ' rem_N    : ', rem_N(layer) 
     644               WRITE(numout,*) 
     645               WRITE(numout,*) ' zd_aon    : ', zd_aon 
     646               WRITE(numout,*) ' zd_eon    : ', zd_eon 
     647               WRITE(numout,*) 
     648               WRITE(numout,*) ' AoN       : ', cbu_i_bio(jn_aon,layer) 
     649               WRITE(numout,*) ' eoN       : ', cbu_i_bio(jn_eon,layer) 
     650               WRITE(numout,*) 
     651               WRITE(numout,*) ' N_C_alg   : ', N_C_alg(layer) 
     652               WRITE(numout,*) ' N_C_det   : ', N_C_det(layer) 
     653            ENDIF ! ln_decoupNC 
    394654 
    395655         ENDIF ! ln_write_bio 
    396656 
    397657      END DO ! layer 
     658    
     659 
    398660 
    399661      !------------------------------------------------------------------------ 
     
    408670     &                  .AND. flag_active(jn_alk) ) THEN 
    409671 
    410             zd_dic1   = - zncp * ddtb 
     672            zd_dic1   = - zncp * ddtb                  !! Bizarre 
    411673            zd_dic2   = - cbu_i_bio(jn_dic,layer) 
    412674            zd_dic    = MAX( zd_dic1 , zd_dic2 ) 
    413675 
    414             zd_alk1   =   no3_c * zncp * ddtb 
     676            zd_alk1   =   no3_c * zncp * ddtb           !! Bizarre 
    415677            zd_alk2   = - cbu_i_bio(jn_alk,layer) 
    416678            zd_alk    = MAX( zd_alk1 , zd_alk2 ) 
Note: See TracChangeset for help on using the changeset viewer.